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Static pressure distribution in the free turbulent jet 


By DAVID R. MILLER* and EXNAWARD W. COMINGS 


Purdue University, Lafayette, Indiana 
(Received 5 March 1957) 


SUMMARY 
Measurements of mean velocity, turbulent stress and static 
pressure were made in the mixing region of a jet of air issuing 
from a slot nozzle into still air. ‘The velocity was low and the 
two-dimensional flow was effectively incompressible. ‘The results 
are examined in terms of the unsimplified equations of fluid 
motion, and comparisons are drawn with the common assumptions 
and simplifications of free jet theory. Appreciable deviations 
from isobaric conditions exist and the deviations are closely 
related to the local turbulent stresses. Negative static pressures 
were encountered everywhere in the mixing field except in the 
potential wedge region immediately adjacent to the nozzle. 
Lateral profiles of mean longitudinal velocity conformed closely 
to an error curve at all stations further than 7 slot widths from 
the nozzle mouth. An asymptotic approach to complete self- 

preservation of the flow was observed. 


1. INTRODUCTION 

Few turbulent flows have received more experimental and theoretical 
study than the simple free jet. Nevertheless, very little appears to be 
known of the role of the mean static pressure in the jet flow. It is commonly 
assumed (Alexander ef al. 1953; Liepmann & Laufer 1947; Pai 1954; 
Schlichting 1955) that the mean static pressure gradient in the direction 
of the jet flow is everywhere negligibly small compared to other mean forces 
on the fluid. Although this assumption seems to be justified in the approxi- 
mate calculation of mean velocity distributions, a consideration of the depar- 
tures from it may lead to a better understanding of the structure of the flow. 

Pitot tube static-tap pressure measurements have been reported for the 
coaxial water jet by Viktorin (1941), for the round compressible air jet 
by Warren (1955), and for incompressible air flows in ducts by Fage (1936) 
and in two-dimensional wakes by Fage and by Schlichting (1930). ‘These 
measurements, while revealing significant departures from isobaric flow, 
were too fragmentary for detaii- interpretation in terms of the basic 
equations of fluid motion. 

In this investigation, the static pressure, mean velocity and turbulent 
stress fields of a two-dimensional, incompressible, free turbulent air jet 
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were measured and analysed. ‘The selection of a two-dimensional flow 
permitted the use of a static pressure probe insensitive to local mean 
flow direction and turbulent velocity fluctuations in the plane of mean flow. 
Although mean velocity data have been published for the two-dimensional 
jet by Forthmann (1934) and by Reichardt (1941), the static pressure and 
turbulence measurements reported here are believed to be unique. 


2. "THEORY 
‘The geometry of the two-dimensional jet is indicated in figure 1. All 
solid boundaries, shown in cross-section, are sufficiently extended in the 
perpendicular z-direction to ensure two-dimensional flow throughout the 
region of measurement in the plane z = 0. Geometrical symmetry of the 
boundaries requires functional symmetry of all scalar flow variables about 


the centreplane (y = 0). 











y, 


— — MOMVOWhS® 


EXPERIMENTAL 
NOZZLE 





Figure 1. Nozzle and flow geometry. 


Equations of motion 
Under the restrictions placed on the flow (two-dimensional, 
incompressible, steady on the average), the appropriate equations of motion, 


with the customary notation, are: 


ee ou ov 
(Continuity) a i 0, (1) 
(x-Reynolds) x (pU? + pu’? + p)+ = (pUV —r+pl) = 0, (2) 
a a. eee 
(y-Reynolds) 5 (pV?+p0%+p)+ 5 (pUV-7-nh) =0, (3) 
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where 7 = —pu’'v’ is the turbulent shear stress and € = dV/ 0x—dU/éy is the 
mean vorticity. Equations (2) and (3) are easily obtained from more 
common forms of the Reynolds equations given by Goldstein (1938) and 
by ‘Townsend (1956) with the aid of (1). 


Boundary conditions 

Those variables which are odd functions of y(V, UV, 7, ¢) are zero at 
y = 0); variables which are even functions in y (U, U?, V?, u, v2, p), have 
zero derivatives in the y-direction at y = 0. In addition, all variables and 
their directional derivatives of all orders approach zero as |y| approaches 
infinity. (In the actual experiment there was a low velocity recirculation 
of air from the downstream to the upstream portions of the jet and the region 
of flow was finite in the transverse direction.) It is convenient to define two 
additional boundaries, y = +A(x), outside which the flow is essentially 
nonturbulent. On and outside these boundaries the turbulent stresses 


(pu'*, pv", 7) and their derivatives are zero. 


Integral forms of the equations 
Making use of the boundary conditions, integration of (1), (2) and (3) 
in the y-direction gives: 


r OU : 

| ~— dy+V =0, (4) 
JQ OX 

| = (pU? + pu’? tp) dy+pUV —7+pE = 0, (5) 
/ 9 % 

[ = (pUV—1—ph) dy +pV*+ pe) +3-p.=0. (6) 


In this investigation, dimensionless forms of these three equations were used 
to calculate the distributions of V, 7 and v” in the fully turbulent region from 
measured distributions of U, u’? and p. 
Taking infinity as the upper limit in (4) and (5) gives the total mass and 
momentum flux integrals ; 
| Udy = I (constant), (7) 
“0 


| (pU?+ pu’? +p) dy = (constant). (8) 
“ @ 
Simplifications used in the theory of fully turbulent flow 
Equations (1), (2) and (3) contain six dependent variables, and solutions 
of the problem cannot, therefore, be obtained by wholly analytic means. 
A number of simplifications based on experiment, assumption or hypothesis 
have appeared in the literature, and are described in the remaining part of 
this section., Heretofore, the lack of sufficiently complete empirical 
information has prevented a critical evaluation of some of these 
simplifications, particularly those involving static pressure assumptions. 
A2 
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1. Reynolds number similarity. \t has been observed by Townsend (1956) 
that in a fully turbulent flow there exists a region including almost all of the 
How in which the direct action of viscosity on the mean flow is negligible, 
i.e. in which 

pUV—7z| > |\u¢|_ for |y A(x). (9) 


The regions of laminar flow outside the turbulent jet are excluded. 


2. Self-preservation. Mean velocity measurements in the fully turbulent 
region of the free jet by Reichardt (1941) reveal that transverse distributions 
of U retain the same functional form at all downstream locations while 
changing only in scale, i.e. 

U = U.f,(n), where n = y/d. (10) 
The scale factors U’. (mean velocity on the centreplane) and 4 (jet width) 
are functions of x alone. In this investigation 6 is defined by the relation 


b(x) = {. (U/U,)° dy. (11) 


‘Yownsend (1956) discusses the hypothetical flow wherein all flow variables 
of equations (1), (2) and (3) (except the vorticity) are self-preserving in the 


same sense as U, 1.e. 





g* = UF ffm), etc. (12) 
He shows that complete self-preservation is theoretically possible in the 
two-dimensional free jet and that necessary conditions are: 
db/dx = c(constant), or b = c(x— Xp), (13) 
and U? ~ (x—xy)-*, (14) 
where kis aconstant. Mean velocity measurements in the two-dimensional 
air jet by Forthmann (1934) and Reichardt (1942), show excellent agreement 
with these conditions and give values for the constants, 
c='0075 and k = 1-00. (15) 
However, measurements of turbulent stresses in the round jet by Corrsin 
(1943) do not exhibit self-preservation within 40 nozzle diameters of the 
nozzle although a tendency toward self-preservation is apparent. ‘This 
has led ‘Vfownsend to conclude that complete self-preservation may be an 
asymptotic condition not valid over the usual range of measurement. 
It should be noted that the nonturbulent flow surrounding an actual jet 


cannot be fully self-preserving. 


3. x-Reynolds equation. In calculations of mean velocity distributions, 
it is commonly assumed, as by Pai (1954) and Schlichting (1955), that 
gradients in the x-direction of the turbulent stress pu’? and the mean static 
pressure are negligible compared to the other forces acting on the mean flow, 
1.e. that 





, (16) 











‘yn 
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When, in addition, the viscous stress is neglected, this assumption leads to 
simplified forms of the x-Reynolds equation, 


5 (pU*)+ 


and of the momentum integral (8), 


(pUV - 72) | 2 (17) 


4 
Ox 


| pU? dy = pbU; =a ip (18) 
0 
4. Hypotheticalindependent equations. Self-preserving analytic solutions 
for U, V and7zcan now be obtained when an independent relation between 
one or more of these variables is assumed. ‘The new relation must be 
independent of equations (1) and (17), which are also used in the calculation. 
A number of such solutions (Gortler 1942; Reichardt 1941; ‘Tollmein 
1926) have been based on such independent relations as the ‘ mixing length’ 
hypotheses of Prandtl (1942), ‘Taylor (1938), and von Karman (1934). 
‘I'wo are presented here for later comparison. 


(a) Gortler’s exchange coefficient 


Arr 








Hypothesis : T = p7(x) ss 

Solution for U: U = U, sech*(2n). (19) 
(6) Reichardt’s thermal conduction analogy 

Hypothesis : oe = A(x) — 

Solution for U: Un expt — hay"). (20) 


Several authors (Corrsin 1943, Liepmann & Laufer 1947, and ‘Townsend 
1956) have noted the shortcomings of these simple hypotheses. 


5. y-Reynolds equation. ‘'Yollmien (1926) assumed that the turbulent 
y-stress (pv) is negligible in the y-Reynolds equation, so that 


A A 


— (pV2+5)4 x (pUV —7) = 0), (21) 


and then solved for the centreplane static pressure ; 

P. = 0-263c?pU-. (22) 
Values of U, V and 7 needed in the calculation were obtained in the manner 
of the preceding paragraph using Prandtl’s first mixing length hypothesis. 
Since c = 0-075 in the two-dimensional air jet, he concluded that the static 
pressure can be neglected in both Reynolds equations. ‘This conclusion 
is in direct contradiction to a theoretical result described by ‘Townsend 
(1956). A dimensional analysis of the magnitudes of the various terms in 
the y-Reynolds equation led ‘Townsend to the conclusion that the lateral 
gradients of the turbulent y-stress and static pressure are the dominant 
terms and that the integrated form of the equation may be written 

p + pv” = p,(x), (23) 


where p,(x) is a function of x and equals the static pressure ‘ outside the flow’ 
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3. EQUIPMENT AND PROCEDURE 

Flow system 

The flow system is illustrated in figure 2. The nozzle, of aspect ratio 
40:1 (20in. by 0-51n.), was specially designed for two-dimensional 
investigations. Control of the turbulence level was provided by two 
15 mesh screens and a honeycomb, followed by a two-stage 56: 1 contraction 
in area. ‘The jet emerging from the nozzle was allowed to mix laterally 
with still air. ‘Ceiling’ and ‘floor’ boundaries helped to maintain the 
two-dimensional character of the jet. All measurements were taken in 
a plane parallel to and midway between them. The floor consisted of 
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Figure 2. Flow system. 


removable plywood sections which provided access for the probes. A 
lead-screw traversing mechanism, located beneath the floor, supported and 
positioned the probes. It permitted manual positioning of the probes 
in the (x, y)-plane with an estimated accuracy of 0-002 in. ‘The centrifugal 
blower was driven at 750 rpm (except during hot-wire calibration) by a 
three-phase induction motor. All measurements were made at a nozzle 
mouth velocity U’,. of 72 ft./sec, corresponding to an exit Reynolds number 
(aU,/v) of 1-78 « 10%. 


Static pressure equipment 


Figure 3 shows the design of the static pressure probe. Reflections of 
the nozzle mouth on the polished surface of the disk were used in aligning 
the probe accurately in the plane of the flow. All pressure readings were 
made relative to a reference tap exposed to stagnant atmosphere far from 
the mixing region. A Prandtl-type micromanometer (Flow Corporation 
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Model MM-2) was used in all pressure measurements. This instrument 
has a range of 0-2 in. and an accuracy of + 0-0004 in. of manometer liquid 
(n-butyl alcohol). No corrections (Fage 1936; Goldstein 1936) were 
applied to the pressure measurements to account for the effect of the z-velocity 
fluctuations on the probe pressure. The maximum error in the measure- 
ments is estimated to be of the order of 5%. This estimate is based on the 
assumption of equality of the intensities of the different fluctuation velocity 
components on the centre-plane of the jet. Equality has been observed on the 
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Figure 3. Static pressure probe. 


centre-line of a round jet by Corrsin (1943) and on the centre-plane of a two- 
dimensional wake by Fage (1936). ‘The values of v® shown in figure 9 
were calculated from measurements of u’2, U and p and were found to 
approach wu on the centre-plane at large distances from the nozzle, as 
would be expected in agreement with the observations of Corrsin and of 
Fage. ‘The measured / is the most important term in this calculation and 
significant errors in the measurement of the static pressure would be reflected 
in corresponding deviations from the equality of uv and v”. 
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Hot-wire equipment 

‘I'he hot-wire anemometer constructed for this study was of the constant 
temperature type and was very similar in design and performance to the 
anemometer described by Laurence & Landes (1952). A Rubicon potentio- 
meter and an electronic mean-square circuit were used in the measurement 
of the DC and AC components of the wire heating current. ‘The tungsten 
hot-wire used to measure both U and wv” had a diameter of 0-0002 in. and an 
active length of 0-08 in., and was supported in a vertical position in the flow. 
Mean velocity hot-wire calibrations were carried out before and after each 
series of hot-wire measurements. An impact tube, mounted in a probe 
with static taps in opposing parallel walls, was used as the standard of velocity. 
Simultaneous readings of the velocity head and the wire heating current 
at two or more velocities constituted a calibration. All calibrations were 
made in the low-turbulence region at the nozzle mouth. ‘The limiting 
factor in the accuracy of measurements was the observer’s ability to judge 
the mean value of a randomly fluctuating visual signal. It is estimated that 
hot-wire mean velocity measurements were subject to a maximum error of 


+3% of U,. 


4. RESULTS AND DISCUSSION 
Hot-wire and static pressure traverses, made at 11 lateral stations and 
along the centre-plane, covered the entire region of turbulent flow out 


a distance 40a (20 in.) from the nozzle. ‘lhe measurements, consisting 


+ 


uC 


Figure 4. Jet width. 


usable dimensionless quantities by a digital computer (Electrodata Model 
203). All lateral profiles were ‘folded’ before plotting by averaging 
computed values on each side of and equidistant from the centre-plane. 


of voltage, current, micrometer and scale readings, were reduced to the 


Only minor asymmetries were apparent in the full profiles. 

Figures 4 through 10 present the main experimental results in dimension- 
less form. Complete tables of data for the single jet reported here and for 
a double jet system are available elsewhere (Miller 1957). ‘The jet width 
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is shown as a function of downstream position in figure 4. The width is 
approximately constant for the first four slot widths from the nozzle. 
Farther downstream, a transition to a linear spread is observed. From 
x/a = 7 to the farthest measurement station, the jet width is well described 
by equation (13), least-mean-square values of the constants being c = 0-0723 
and x,/a i572. 

‘The decay curve for the centreplane velocity (squared) appears in 
figure 5. ‘lhe straight portion of the curve for x/a> 7 hasa slope of — 1-028 
as compared to a slope of —1 predicted by the approximate equation (18). 
The behaviour of 6 and U, with distance from the nozzle defines two 
distinctive flow regimes, the transition flow region extending from the nozzle 
to x/a = 7, and the fully turbulent region lying farther downstream. 
‘These two regions are considered separately. 
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Figure 5. Centreplane velocity decay. 


Fully turbulent region 

Figure 6 is a composite plot of the U-profiles measured at x/a = 10, 20, 30 
and 40. In this and other figures pertaining to the fully turbulent region, 
U.and 6, both functions of x, are used to normalize the data. All experi- 
mental points lie on the single curve, within the limit of error, indicating a 
high degree of self-preservation of the mean velocity profiles. ‘The plotted 
points in the low velocity region do not represent U truly since the single 
hot-wire does not distinguish between U and V. ‘The points actually 
represent ,/(U*+17)/U.. Making use of this fact and equation (4), the 
true values of U can be computed. Both the corrected and uncorrected 
profiles are shown in figure 6, as well as the square of the corrected profile. 
The smoothed and corrected experimental curve is seen to be in good 
agreement with points on the error curve of Reichardt’s analysis (see (20)). 
Gortler’s profile (see (19)) deviates widely from the data for values of /b 
greater than 1-8, indicating that the eddy viscosity is not constant across 
a section but decreases with distance from the centre-plane. 
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Turbulent x-stress profiles appear in figure 7. ‘The most striking feature 
of this plot is the merging of the curves into a single envelope curve in the 
yuter region of mixing. Partial self-preservation is thus exhibited with 
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Figure 6. Mean velocity profiles 


a detinite tendency toward complete self-preservation as the distance from 
the nozzle increases. Each curve of figure 7 is well correlated by a 
summation of two error curves. From the trend in the error curve 
parameters with distance, it is reasonable to expect that the asymptotic 
profile is approximated by the empirical relation 

(u’?/U?), _ ,, == 0-17 exp( — 0-517?) — 0-10 exp( — 1-157). (24) 
The outer boundary of turbulence is seen to fall at y/b = 4, i.e. A = 46 
in the fully turbulent region. 

Static pressure data are plotted in figure 8. ‘The gradual transition of 
the profile with longitudinal distance is apparent. Partial self-preservation 
is not observed in these curves. However, the trend in the shape of the 
curves is toward a stable functional form, varying only in amplitude as the 
distance from the nozzle increases. ‘The two farthest profiles are well 
described by the simple error curves 

[(P—P,)/pU Fra = 30 = — 9°054 exp(— 0-30n?), (25 a) 
[(p —p:)/pU 7] pia — 49 == - — 9-060 exp( — 0-30). (25 b) 
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Figure 8. Static pressure profiles—fully turbulent region. 


The correspondence is better for the profile at «/a = 40, indicating that 
functional self-preservation is achieved close to this station. Note that p; 
is neither zero nor constant with x as is sometimes assumed. Because it is 
not self-preserving, it has been subtracted from p in equation (25). 








12 David R. Miller and Edward W. Comings 


A comparison of the empirical coefficients of equations (25a, b) with 
the corresponding centre-plane values of u’*/U? of figure 7 reveals that the 
quantities are very nearly equal. If it is assumed that this equality persists 
for larger values of x, the asymptotic solution for the static pressure is 

[((p—p,)/pUF]e = 0 = — 9-07 exp( —0-30y?), (26) 
where the coefficient is obtained from equation (24). 

‘Tollmien’s solution for the static pressure on the centre-plane (see (22)) 
gives p,/pU? = 0-00137 when c = 0-0723, whereas the measured values are: 
v/a = 10 20 30 40 

pip? = —0-0259 —0-0389  —0-0498 — —0-0582. 


The lack of equality indicates that Tollmien’s basic assumption (21) is 


grossly in error. 
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Figure 9. Turbulent v-stress and shearing stress profiles—fully turbulent 
region. 


Using the data of figures 4-8 in equations (4), (5) and (6), values of 
V,«’2and 7 were computed. ‘The latter two variables are plotted in figure 9 
together with p and uw’ for comparison. ‘There is a close negative corres- 
pondence, improving with v, between Pp and pu? 1.€,. 


79 


p+ po? = B(x). (27) 
This relation is identical to the one put forward by ‘Townsend (see (23)), 
and it supplies experimental verification of the validity of his neglect of 
other terms in the y-Reynoldsequation. All remarks made above concerning 
self-preservation of the static pressure profiles apply also to the turbulent 
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y-stress profiles. With a change in sign, equations (25 a, b) are equally good 
descriptions of the v’?/U? profiles. Figure 9 reveals clearly the tendency 
for v2 to equai uw” on the jet centreplane. 

Self-preservation of the turbulent shear stress profiles is seen to exist 
at x/a = 20, 30 and 40. ‘This is to be expected since the other significant 
terms in equation (5), U? and UV, are both fully self-preserving. ‘The 
fact that these profiles do not approach zero at y/b = 41s taken as an indication 
of errors of measurement and calculation. ‘lhe magnitude of the error 
at the limit of the region of turbulence, where 7 should approach zero, is 
particularly sensitive to the value of & used in the calculation. A value 
closer to 1-000 than the measured value 1-028 would result in a reduction 
of the error. 

In order to test the assumption underlying equation (16), the values of 
the two sides of the equation were calculated at a number of different points 
in the flow and are set out in table 1. ‘The implication of the assumption 
is that the numbers appearing in the last column of the table should be much 
less in absolute magnitude than the numbers in the adjoining column. 
This is seen to be the case at most points; however, there exist isolated 
regions where the y-derivative passes through zero and the assumption is 
invalid. Even in these regions equation (17) is not seriously in error, 























vla | yla — © (pUV—r+pl) | —, 2 (pw’?+3) | 
pl 2 dy | pl 2 Ox | 
10 0 0-048 | 0-0002 | 
10 1 0-010 | 0-0013 | 
10 | 2 0-001 | 0-0014 | 
10 3 (0-000 | 0-0001 | 
20 0 0-012 0-0003 

20 1 0-003 | (0-007 
20 2 0-004 0-0005 
20 3 | 0-002 0-0007 
30 0 | 0-007 | 0-0001 | 
30 El 0-004 | 0-0003 | 
30 ze | 0-000 0-0004 
} 30 | 3 | 0-002 | 0-0002 | 
| 40 0 | 0-003 | 90-0000 | 
40) | 1 | 0-003 | 0-0001 | 
Aye 2 | 0-000 0-0002 
40 | 3 | 0-001 | 0-0002 

| | 





Table 1. Test of approximation underlying equation (13). 
\ ] ) 


because of the small magnitude of the neglected terms. ‘This is a fortunate 
situation, due largely to the fact that pu’ and p are of the same order of 
magnitude and of opposite sign. ‘Their sum and its derivatives are con- 
sequently much smaller than either term and its derivatives alone. ‘The 
x-Reynolds equation simplification leading to (16) is thus justified in free 


jet flow. 
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Transition region 

Lateral distributions of U, u’2 and p were measured at x/a = 0, 0-6, 
1,2,3,4and6. ‘Thecomputation procedure for obtaining the corresponding 
distributions of V, v’ and 7 was not used because of the relatively large errors 
involved in calculating derivatives in this region of large gradients. 

The region of low turbulence level (+/(u’?)/U < 0-05) immediately 
downstream from the nozzle was contained in a wedge whose vertex fell 
on the centre-plane at «/a = 2:5 and whose base coincided with the nozzle 
mouth. A positive static pressure ridge, straddling the centre-plane and 
decaying rapidly with x, was found within the ‘potential wedge’. ‘This 
ridge is attributed to the persistence of the pressure distribution created 


within the nozzle. 




















Figure 10. Turbulent x-stress and static pressure profiles—transition region. 


Outside the potential wedge, all measured static pressures were negative. 
Particularly striking were the static pressure trenches appearing on each 
side of the potential wedge starting at x/a = 0-8 and finally merging at the 
centreplane at x/a=5. Static pressure measurements reported by 
Viktorin (1941) and Warren (1955) reveal similar side trenches in the coaxial 
water jet and round compressible air jet, respectively. Several of the static 
pressure profiles are plotted in figure 10, the central ridge and side trenches 
being clearly visible. 

The corresponding turbulent x-stress profiles of figure 10 show the 
development of turbulence in the high shear regions on either side of the 


vi 


Static pressure distribution in the free turbulent jet 1 


potential wedge. Except in the potential wedge, a close negative corre- 
spondence is noted between the static pressure and turbulent x-stress profiles 
at x/a=1 and 2. ‘This correspondence is seen to deteriorate at stations 
farther downstream. 


5, CONCLUSIONS 

A negative static pressure field exists throughout the turbulent region 
of the two-dimensional free jet. Its magnitude and directional derivatives 
are of the same order of magnitude as those of the x and y components of the 
turbulent stress and the turbulent shear stress. 

Static pressure is of minor importance in the x-Reynolds equation of 
motion for two reasons; (a) the longitudinal mean momentum flux pU? is 
on the order of 20 times the magnitude of the static pressure, and (b) the 
longitudinal turbulent stress pu’? is of the same order of magnitude as the 
static pressure but of opposite sign, thereby tending to cancel its influence. 
The approximation underlying equation (16) is justified despite the fact 
that it is not correct in certain regions of the flow. 

The static pressure p and the lateral turbulent stress pv’? are the two 
dominant stresses in the equation of lateral motion in the fully turbulent 
region. ‘lollmien’s (1926) analysis and conclusions are therefore invalid, 
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whereas 'l‘ownsend’s (1956) analysis is confirmed. 

A close correspondence between the static pressure deficiency and 
turbulent stress components (pw near the nozzle, pv in the fully turbulent 
region) leads to the conclusion that where the flow is turbulent the static 
pressure deficiency is a manifestation of the presence of turbulence. 

Self-preservation of the mean velocity and turbulent shear stress profiles 
linear spread of the jet, and negligible viscous shear stress, were all observed 
in the fully turbulent region far from the nozzle, as expected from past work. 
A tendency toward self-preservation, observed in the static pressure and 
turbulent stress profiles, supports ‘Townsend’s conclusion that complete 
self-preservation of the flow may be an asymptotic condition not valid over 
the usual range of observation. ‘The assumption that asymptotic profiles of 
U, v, and p, are simple error curves and that the asymptotic wu” profile is a 
compound error curve is not inconsistent with the reported data. 

Of the various semi-empirical mean velocity solutions, Reichardt’s 
simple error curve describes the data best. ‘The eddy viscosity (or exchange 
coefficient) is not, therefore, a function of x alone and Gortler’s analysis 
yields a poor representation of the mean-velocity profiles. 
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SUMMARY 


In order to answer some of Proudman’s questions (1956) 
concerning shear layers in rotating fluids, a study is made of the 
How between two coaxial rotating discs, each having an arbitrary 
small angular velocity superposed on a finite constant angular 
velocity. It is found that, if the perturbation velocity is a smooth 
function of 7, the distance from the axis, then the angular velocity 
of the main body of fluid is determined by balancing the outflow 
from the boundary layer on one disc with the inflow to the boundary 
layer on the other at the same value of r. At a discontinuity in 
the angular velocity of either disc a shear layer parallel to the axis 
occurs. If the angular velocity of the main body of the fluid is 
continuous, according to the theory given below the purpose of 
this shear layer is solely to transfer fluid from the boundary layer 
on one disc to the boundary layer of the other. It has a thickness 
O(v'3), where v is the kinematic viscosity, and in it the induced 
angular velocity is O(v"®) of the perturbation angular velocity of 
the discs. On the other hand, if the angular velocity of the main 
body of fluid is discontinuous, according to the theory given 
below the thickness of the shear layer is O(v'*). A secondary 
circulation is also set up in which fluid drifts parallel to the axis 
in this shear layer and is returned in an inner shear layer of 
thickness O(v'?). 

The theory is also applied to the motion of fluid inside a closed 
circular cylinder of finite length rotating about its axis almost as if 
solid. 


1. INTRODUCTION 


Recently Proudman (1956) has investigated the motion of an in- 
compressible viscous fluid confined between two concentric spheres which 
rotate about a common diameter / with angular velocities (2 and 0(1 +), 
where € is small. He found that when the kinematic viscosity v is small the 
circular cylinder C circumscribing the inner sphere and having its generators 
parallel to / separates out two regions with different properties. Outside C 
the fluid rotates as if solid with the angular velocity of the outer sphere while 
inside C the angular velocity is a function of distance from / only, deter- 
mined by balancing the fluid entering the boundary layer on the faster 
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moving sphere with the fluid leaving the boundar, layer of the slower 
moving sphere. According to this solution there must be a return of fluid 
and Proudman inferred that it could occur only near C, but he was unable 
to decide whether the appropriate shear layer would have a thickness 
O(v18), O(v'*), or even more, and did not determine the behaviour of the 
fluid in any of these possibilities. In view of the difficulty of determining 
even the principal properties of rotating fluids in certain problems it is 
desirable that Proudman’s argument be completed so as to eliminate the 
possibility that the flow between the spheres be of an entirely different 
character. 

The problem considered by Proudman is unfortunately too difficult to 
solve completely at present, and so in this paper a simpler problem, namely 
the flow between two coaxial planes rotating almost as if solid, is considered. 
A complete formal solution may be obtained which supports his inferential 
argument when v is small. In addition, if there is a discontinuity at a 
distance a from the axis in the angular velocity of either discs, a shear layer 
of the character required by Proudman is formed near the coaxial cylinder 
of radius a. ‘In this layer any necessary change in the angular velocity of 
the main body of the fluid occurs in a distance O(v!4) while any necessary 
transfer of fluid from one disc to the other occurs in a distance O(v'*). 
Associated with the change in angular velocity there is aiso a circulation 
in which fluid is carried parallel to the axis in a layer of thickness O(v""*) 
and returned in an inner layer of thickness O(v!*), A like problem also 
discussed here is that of two coaxial parallel discs rotating with the same 
angular velocity but enclosed in a coaxial cylinder of radius a rotating 
with a slightly different angular velocity. The boundary layer on the 
cylinder also consists of two parts and on the discs the boundary layer is 
as usual of thickness O(v'?) but only penetrates a distance O(v'*) from 
the cylinder. 


2. Flow between coaxial rotating discs 

Consider two discs in the planes s = +d and zs = —d rotating about 
the = axis with angular velocities 

Q4 QOctFy(r)+Fi(r)t and Q+Qe{F,(r)— Fyr)}, 

respectively, where ar and az denote distances from an along the axis of 
rotation, a is a characteristic length and ¢€ is small. Let the components of 
the velocity of the fluid be (wu, v, w) relative to fixed cylindrical polar axes 
with coordinates (r, 0, ). Then since by symmetry all dynamical variables 
are independent of the azimuthal angle @ we may use the equation of con- 
tinuity to define a stream function #& by 





eaQ) Ous eaQ) ous 
u= —-—=z~, w= =. (2.1) 
r 02 r or 
Further write rv = aQr*? + eaQy. (2.2) 


If « = 0 the solution of the Navier-Stokes equations which satisfies the 
boundary conditions is y=%=0. Hence since € is small it is reasonable 
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to retain only the first order terms in «, when as Proudman has shown the 
Navier-Stokes equations reduce to 


a “2 a “9 9 
a ( X CG” 1 ( 9 ie r4 
SRS chee hee eo (2.3) 
C3 or” Y of On” 
5 , Ob } 02 l C fae 5 
teen eer a Te (2.4) 
Ox 1 op ¥ o} C3 


where R = a@?Q/v. ‘The boundary conditions are that 


%=des/os =0 ats = +d, | 
x=rP{F(r)+F(r)} atz=d, } (2.5) 
¥ = PiF(r)—F(r)} ats = —d. | 


The main interest of the problem stated above lies in its exemplifi- 
cation of the flow arising from more general rotating bodies, such as the 
spheres considered by Proudman (1956). For rotating discs the particular 
case when /’,, /, are piecewise constant is of greatest interest. It is con- 
venient to treat the problem in two parts by setting /'; and F, successively 


equal to zero. 


3. "THE ANTISYMMETRICAL PROBLEM 
Write 


x=r| A(k)J,(kr)sinh «z dk, 





w=r| B(k)J,(kr)cosh az dk. 
‘These functions are solutions of (2.3) and (2.4) if 
BRo? = —(#?—k?)A, (x? —k?)8 = —4R?* o?. (3.2) 
‘There are three roots %,, %, %3, each with positive real parts, of the 
sextic in « which lead to independent solutions. Hence the general 
solution of (2.3) and (2.4) of the form (3.1) is 
x=r> | A,fk)J{(Ar)sinh «, 3 dk, 
8 1/0 a ax 
; 2 2B ; ee) 
"2 gy? _ R2 ( 


b=? s | = A,(k)J,(Rkr)cosh x, = dk. 


K 











> A,sinha,d = fk), > A,(a2—K)sinha,d = 0, | 
s— | 8=1 e 
3 — (3.4) 
dl 42 — k* f 
A, — cosha,d = 0, 
2 ae: OSN &, ¢ 
where rF(r) =r | fo(k)J (rk) dk. 


From the three equations A, can be determined in terms of f;, and a 


complete formal solution found. Of particular interest is the form of the 
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solution when R is large, and during the remainder of the paper we shall 
concentrate on this special case. ‘Then so long as k < R'* the relevant 
roots of the sextic (3.2) are 
oh an + O(k'R*%), to = RV2(1+72)+ O(F?R?), 
a3 = RN2(1—71)+O(R? RU), (3.5) 
It is noted that x and x, are associated with regions of rapid change in the 
axial direction of thickness O(R'!'?) while x, is associated with a region of 
rapid change in the radial direction of thickness O(R'*). It will appear 
later that all three can combine to give another region of rapid change in the 
radial direction of thickness O(R''*). ‘The assumption k < R'? is justified 
only if the integrands in (3.3) are negligibly small when it is not satisfied. 
It will appear below that this is true except within a distance O(R-'?) of any 
irregularity in the angular velocity of either disc. 
Further, 
Rf + 
A, = DRI cosh(Bd/2R) +Rsinh(Rd/2R)’ o 
an kid 2 
A, exp{R'?d(1+7)} = Agexp{R!d(1—72)} = — A, cosh IR’ (3.7) 








provided that, in addition to k < R'*,d> R*. 

If F,(r) isa smooth function tending to zero as r—> «, k?d may be neglec- 
ted in comparison with R in the determination of x and ys. We find that 
when |z|—d = O(R'?) the contributions from A, and Ay, are 


except z 
exponentially small and 
sr 
x FR 
O2R/ Rls ; 
| (She ) valk) dk = 1R"*r*f,(r). (3.9) 





RS (k)J (Rr) dk 3.8 
kif (3.8) 


while oe 
b= oR}, IRE 
In the neighbourhood of z = +d there are, in addition to the contribution 
from A,, also important contributions trom A, and A, leading to boundary 
layers of thickness O(R-?) in which 
y= PP ire *cosx, | 
b = {R-V2,2F,(r)[1 —12 e-* cos(v—j7)], (3.10) 
« = (d—|2\|)R"?. 

‘These results agree with Proudman’s arrived at inferentially. ‘lhe 
outflow from the boundary layer of the slower moving disc must be the 
same as the inflow to the boundary layer of the faster moving disc. Using 
Proudman’s argument this could only be achieved if the angular velocity 
in the main body of the fluid were the mean of the angular velocities of the 


discs at the same value of r, which is in agreement with (3.8). 
If d = O(R-?) a special treatment is required but there is no difficulty. 
Physically it means that the boundary layers on the discs have joined up. 
If F, has a discontinuity, at r = 1 for example, ¢# also has a discontinuity 
and the integral for y fails to converge there. A closer examination of the 
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flow near r = 1 is necessary and for this purpose let us take 


rwF,=-1 (r<1),) (3.11) 


Although the same approach as before may be used with success the easiest 
method is to assume that modifications to the solution already obtained 
are necessary only in the neighbourhood of r = 1. It is then legitimate to 
neglect the operator (1/r) 0/er in comparison with 0?/dr? and _0?/0s* in (2.3) 
and (2.4), and to extend the range of r from (0, «) to(—@, #). ‘The 
problem is now reduced to a form amenable to the sine transformation. 
Write 
ee 


¥ = p. | A,(k)sin k(y — 1)sinha, z dk, 











| . i (3.12) 
2. fie kee | ; 
y= — IR z |. = : A(k)sink(r—1)cosha,z dk, | 
where «, is defined in (3.2), (3.5), and A, in (3.4), except that now 
folk) = 2/ak. (3.13) 
Hence, except when |z|—d = O(R-*?), 
7 sen(r — l) g J [" e’kir-l’ sinh(k33/2R) dk (3.14) 


}_ ,, 2R'* cosh(k?d/2R) + k sinh(k?d 2R){’ 


it being assumed that the integrand is negligibly small when k = O(R!?). 
The integral may be evaluated by contour integration on noting that 
since R is large the poles of the integrand occur at cosh (k?d/2R) = 0, and 





we have 
: 2 sgn(r— 1) iP sin{(22 + 1)z3/2d} _ 
. 372 SqISRIO Ss, vadae rate 
x [e Pn — 2ePn? cos($V38,, —47)| + O(R13), (3.15) 
where — (On + 1)7R)\ 13 
le eee: aad! 1}. 


Similarly, 
sen(r—1) 2sgn(r—1) (—)" cos{(2n + 1)z2/2d} 
| i a 2n+1 

x [ePn + Zen? cos 4V38,]+ O(R%). (3.16) 
Near 2, = dextra terms must be added from A, and A, leading to boundary 
layers of the type described in (3.10). Thus the thickness of the shear 
layer near r = 1 necessary to permit the return of fluid from one disc to the 
other is O{R*) and in it the fluid velocity is O(R™*). 

The particular form (3.11) was chosen for F, because of the simplicity 
of its sine transform (3.13), but it is to some extent artificial since 
the angular velocities of the discs are singular at r=0 and in any 
case the linearization of the equations of motion is invalid there. 
Nevertheless the solution obtained describes the shear layer appropriate 
to a discontinuity in the angular velocities of the discs at r = 1 because its 
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effect is confined to the immediate vicinity of r = 1, and further, outside it, 
# = }R? spn(r— 1), x =o(R-*), (3.17) 
in agreement with (3.8) and (3.9). 

The solution can also be used to complete the description of the flow 
given in (3.8) and (3.9) in the more realistic case when F, is bounded. Let 
F, = U(r)—m(r) (r < 1),) 
=I(r)+m(r) (r> 1), J 
where /, m are smooth bounded functions of r, and write yp, %» respectively 

for the right hand sides of (3.15), (3.16). Then 

x = mr)x9 + OCR), 

% = R17? 1(r) + 7r?m(r) by + OCR) (3.19) 
is uniformly valid everywhere except in the boundary layers of the discs. 
The effect of the approximations made so far may now be established. 
‘The operator (1/r) ¢/er makes contributions O(R ' 349) to % and O(R-! Fy) 
to y. Further, since all the contributions from the shear layer tend to zero 
exponentially outside it the effect of extending the range of the problem 
from 0 <r < © to —% <r < 18 to add exponentially small terms to 
x andy. Hence in general the solutions given above are the correct leading 
terms when R is large. ‘There is a breakdown, however, if that part of the 
range of integration when k = O(R'”) makes a significant contribution to 
any of the integrals. A simple way of estimating when this is likely to 
happen is to examine the behaviour of the approximate forms when 
k = O(R'?). The integrand when F,(r) is smooth will certainly be 


(3.18) 


negligible in virtue of f,(k). Let us consider (3.14) as an example of the 
situation when F,(r) is not smooth. Writing k = k,+7k,, the modulus of 
the integral when k, > 0 and k = O(R'?) is 

l d—|z| _, tia 2] 
Rt R12 exp, — k, r—1|- Riz k — 3k5k, j 
and is not exponentially small only when both |r—1) and d—'z are simul- 
taneously O(R''?), that is, only in the immediate neighbourhood of the 


discontinuity in the angular velocity of the discs. 

Finally Proudman examined the order of magntidue of the non-linear 
terms in the Navier-Stokes equations and from his results we deduce that 
the theory of the flow in the shear layer is valid except aty = 1, s =din 
the limit e->0, R-> «© sothateR—-0O. ‘The condition is more than sufficient 
to ensure the validity of the theory elsewhere apart from the artificial, but 
convenient, case described in (3.11) which fails at r = 0. 

It is possible for F,(r) to be continuous but to have a discontinuity in 
its first derivative. An example is 


Pedr) =¢7 (3.20) 


where y is a constant. Using the same technique as when F was dis- 
continuous we find that the shear layer is again of thickness O(R''*) and 


that in it y = O(R"”), 4 = O(R**). 
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4. THE SYMMETRICAL PROBLEM 


In this case 
x=PFr) atz= +d, (4.1) 


and we write 


s 


yor [ A (k)J ,(kr)cosh «, z dk | 
, (4.2) 


yr 3. 7© a, —k? f 
b= — IR > |, —— A (k)J,(kr) sinha, z dk, 


and the treatment closely paralleis that for the anti-symmetrical problem, 
the only difference being that cosh and sinh are interchanged throughout. 
If F’, is a smooth function it may be shown that 
x =r F\(r)+O(R?) and = O(R4), (4.3) 
everywhere confirming Proudman’s argument. Near the discs there are 
boundary layers in which the change in y is O(R"!”) and the change in & is 
O(R"). The flow near a discontinuity however introduces novel features. 
For defining F;, as in (3.11), 
sgn(r — 1) [a e'* -l' cosh(k3z/2R) dk 


— 7 5 | | , 2R'* sinh(k?d/2R) + k cosh(k?d/2R) 





(4.4) 


in the immediate neighbourhood of r = 1 except if = —d = O(R-?), when 
extra terms, as in (3.10), must be added. Since R is large the poles of the 
integrand occur at sinh(k#d/2R) = 0 and at k°d+R'2=0. Evaluating 
the residues at the relevant poles we have 


x = sgn(r— 1) E —exp{ — Ri4d-1?\7 —1|} - 














m 
2 2. (—)" cos(n73/d) , 
— 2 ees 9) 2 i 20 1 2 5 
3R1U6q! 5 2, (2nz)? ss ™ . COS (2 \ IP, 3 ff (4 3) 
*here (2Rnz) is 
= here VY, = 2 7 : r—1| ° (4.6) 
¢ 
Similarly 
p= — spiky sgn(r— l)exp{ — R™4d-42\7r — 1}} + 
sen(r—1) 2 (—)" sin(nz2z/d) 8 
bs aa on + De—nl2 eos(1V3y_)]. ' 
3nRIe : le e cos(}v3y,,)]. (4.7) 


In this case therefore the shear layer has a thickness O(d!?R-!*), con- 
firming one of Proudman’s suggestions. A shear layer of this type is 
necessary to allow y to change sign. Asa result fluid drifts from the central 
plane towards the discs with velocity O(R-'*) and the return of fluid from 
the discs takes place in an inner boundary layer of thickness O(R~'*) and 
the azimuthal velocity that is induced is O(R“'®). 

As d-> « the contribution to y in (4.3) from the first two terms tends 
to zero and that from the series is O( R-'®) except in the immediate vicinity 
of z|=d. ‘Thus if a disc is rotating in an unbounded fluid the angular 


velocity of the fluid cannot have a finite discontinuity in the limit of zero 
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viscosity, for the shear layer it would generate when the viscosity is small 
could not be thin. ‘This result is to be contrasted with the related problem 
of the flow due to the slow motion of a body of revolution along the axis of 
an unbounded rotating fluid. In the special case of a circular disc of unit 
radius Morrison & Morgan (1956) have obtained a complete solution. 
They showed that when v is small the effect of viscosity may be neglected 


save near the disc and the cylinder r= 1. In the inviscid solution the 
perturbation of the angular velocity is zero outside and O{(i —r) 1} inside 
the cylinderry = 1. Ashear layer is necessary to smooth out this singularity ; 


it is controlled by the same set of equations as in the problem discussed 
above and its thickness is O(R''*). Paradoxically the shear layer can 
smooth out the infinite discontinuity in Morgan’s problem but cannot 
smooth out the finite discontinuity in the present problem. It would be 
interesting if Morgan’s problem could be solved for a fluid bounded above 
and below the disc to see whether there is any fundamental change in the 
character of the shear layer. 

A parallel treatment, defining /, as in (3.17), shows that, if the angular 
velocity of the discs has a discontinuous derivative, in the shear layer 

= O(R-**) and the increment of x is O(R“*). 


5, ROTATING CYLINDER 
Another example in which one of the friction layers has a thickness 
O(R'*) is provided by the flow inside a circular cylinder which is almost 
rotating like a rigid body. Suppose it has its plane surfaces in the planes 


z= +d, the zs axis for its axis of symmetry and of rotation, and unit radius. 
Let the angular velocity of the plane surface be Q and of the curved surface 
be Q(1+¢e). ‘Then proceeding as before we have for boundary conditions 
A] 5 
OW 
b=~=y=0 ats td (<r 1), | 
O38 | on 
Ous (5-1) 
boa =yxy-l=0 atr=1 (2 d). | 
or | 


It is likely, and will be confirmed later, that the perturbation in the 
velocity components from those in a state of uniform rotation with angular 
velocity Q will be confined to the neighbourhood of r = 1. Accordingly 
we neglect the operator (1/r) 0/ér in the linear differential equations satisfied 
by y and ¢ and extend the range of r to(— «, 1). 

Then write 


yo > S Ae a eee, 2 
ead 
whence 
l :. a” — ke ' 2 
ys = “Re 2 SY A,,, A hn” sin a5, 2; (5.2) 
7 - ' 


Sn 


where «,,, is any one of the cs roots of 


(02 —k2)3 + 4R2x2 = 0 (5.3) 





2 
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with positive real parts, and the summation is over all acceptable values of 
k,. When R is large and k,, < R"? the relevant roots of (5.3) are 

a, = RIZR, R'?(1 +7), a, = R2(1 —2). (5.4) 
‘he acceptable values of ad which must have positive real parts, are 
determined by the conditions at |z| = d. ‘These lead to 


> A,,, cosa,, d = 0, > «;,, Acosa,,, d= 0, 
] s ] 


— A. sineg,, d = 0, 


n Si 





which have a non-zero solution for A,, only when the determinant of the 
coefficients vanishes. For large R this means that either 


sin(k? d/2R) = 0, k, 40, that is, 


/ 


2Rn7\13 2 a\ 1/3 ? a\ 1/3 
k, i n ( Rn ) (4-443!) ot (- in | (4—4v31), (5.6) 











d d ee 7 * 
where 7 is a positive integer, or k, = (R/d?)!4. When this condition is 
satisfied, 
k3 d : 1/2 ‘5 ae ibe 
1, cos 5% = 2A, cos R'*(1+7)d = —2A,cos R'*(1 —2)d. (5.7) 


The contributions from A, and A, may be neglected except in the immediate 
neighbourhood of |z| = d and constitute a boundary layer in which the 
perturbation components in the main body of the fluid are brought to zero 
at the plane faces. Hence elsewhere, substituting the permissible values 
of k, into (5.2), we find that 

= Bexp{ — R"4d-42(1 — 1); 


1 “ cos(n73/d) 








* QiBRU & ~ (Qan)P8 (C,, "+ D,, € 7"? cos(3v3y, +7,)§, (5-8) 
y= — ...# exp{ — RM4d-42(1 —r)5 4 
io : 
| > sin(ums/d r 
¥ 27 R! 2 p2 ’ nN | iC ia é a D, é — cos(5 \ IY n re Yn a 477)}, 
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where y,, is defined in (4.6), and B, C,,, D,, 7, are constants. From the 
boundary conditions at r = 1 where y,, = 0, we have 


B= j, C,,+C,, cosy, = 90, (since x = 1), (5.10) 
Ors 

C,,—D,, cosy, = O(R-™!”), (since ie 0), (5.11) 

C,,+ D,, cos(y, +47) = 2(—)"*1, (since % = 0). (5.12) 


Thus the leading terms are 
: ’ 4 + ef 
x = exp{— RY4d-*"(1 —1)5 4 312q13 RUG 2 
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‘hus again the shear layer near ry = 1 is in two parts. ‘The main part has 
thickness O(R-!*) and in it y is reduced from unity to zero. However to 
do this + and 0xs/ér must be given non-zero yalues on the curved surface of 
the cylinder. An inner boundary layer whose thickness is O(R-%) is 
therefore necessary to bring these to zero on the wall. This in turn leads 
to an inner layer for y but the change in y which results is only O(R®). 
Finally near the planes s = d additional boundary layers like (3.10) occur 
to bring the dynamical variables to zero at |z| = d. It is noted that these 
will extend for a distance O(R-"*) inwards from r = 1 and O(R-'?) from 
z_ = d while the basic assumption of this paper, that k? < R 1s not justified 
if z|—d and 1—,r are both O(R-?!?). 
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SUMMARY 

The instability of small gas bubbles moving uniformly in 
various liquids is investigated experimentally and theoretically. 

The experiments consist of the measurement of the size and 
terminal velocity of bubbles at the threshold of instability in 
various liquids, together with the physical properties of the 
liquids. ‘The results of the experiments indicate the existence 
of a universal stability curve. The nature of this curve strongly 
suggests that there are two separate criteria for predicting the 
onset of instability, namely, a critical Reynolds number (202) 
and a critical Weber number (1:26). The former criterion 
appears to be valid for bubbles moving uniformly in liquids 
containing impurities and in the somewhat more viscous liquids, 
whereas the latter criterion is for bubbles moving in pure, 
relatively inviscid liquids. 

The theoretical analysis is directed towards an investigation 
of the possibility of the interaction of surface tension and hydro- 
dynamic pressure leading to unstable motions of the bubble, 
i.e. the existence of a critical Weber number. Accordingly, the 
theoretical model assumes the form of a general perturbation in 
the shape of a deformable sphere moving with uniform velocity 
in an inviscid, incompressible fluid medium of infinite extent. 
The calculations lead to divergent solutions above a certain 
Weber number, indicating, at least qualitatively, that the interaction 
of surface tension and hydrodynamic pressure can result in 
instabilities of the bubble motion. 

A subsequent investigation of the time-independent equations, 
however, shows the presence of large deformations in shape of 
the bubble prior to the onset of unstable motion, which is not 
compatible with the approximation of perturbing an essentially 
spherical bubble. This deformation and its possible effects on 
the stability calculation are therefore determined by approximate 
methods. From this it is concluded that the deformation of the 
bubble serves to introduce quantitative, but not qualitative, 
changes in the stability calculation. 


* Now at Cornell Aeronautical Laboratory, Buffalo, New York. 
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1. INTRODUCTION 

Very small bubbles rise in a rectilinear path in a body of liquid at rest. 
It is found, however, that if one forms larger bubbles in liquids of relatively 
small viscosity, a critical size for each liquid is reached at which the bubble 
suddenly assumes an oscillatory trajectory. ‘The actual path is most often 
helical, although the bubble is occasionally seen to oscillate in a single plane 
as it rises. ‘lhe bubbles in this region are seen to be slightly deformed ; 
they seem to have the shape of oblate spheroids. Still larger bubbles become 
more deformed and rise rectilinearly, but with a rocking motion about the 
axis in the direction of their motion. Finally, a size is reached at which the 
bubble rises rectilinearly maintaining a spherical shape on its upper surface 
but a very irregular, fluctuating shape on its lower surface. ‘These bubbles 
have been called ‘spherical-cap’ bubbles. For more viscous liquids, such 
as mineral oil, the bubbles simply become more and more deformed as the 
size is increased, with no oscillations in shape. Accordingly, all bubbles 
in such liquids rise rectilinearly. 

It is the purpose of this paper to examine experimentally the conditions 
necessary for the onset of oscillations in relatively inviscid liquids, and to 
develop a theory to predict these ‘critical’ conditions. 

Most of the previous theoretical research has been limited to solutions 
for the drag of rigid and fluid spheres moving slowly and uniformly in an 
infinite medium. Usually, the inertial terms in the Navier-Stokes equations 
have been neglected. Notable exceptions to this trend of research are the 
solutions for the shape and motion of spherical-cap bubbles by Davies & 
‘Taylor (1950), and for the drag of rigid spheres at a Reynolds number of 
about 100 by Kawaguti (1955). Kawaguti studied the stability of the flow 
about spheres and found instability at a Reynolds number of 51, which is 
approximately one-half of the observed value for rigid spheres. His 
treatment is not applicable to fluid spheres because of the slip boundary 
condition required. Recently, Saffman (1956) obtained relations for the 
terminal velocities of spiralling bubbles in water by assuming the flow near 
the front of the bubble to be inviscid and considering the distribution of 
pressure in the vicinity of the stagnation point. By treating the planar 
oscillating bubbles in the same way, he arrived at an equation which 
determines the stability of the rectilinear motion of bubbles in water. ‘The 
value of the critical Weber number U(pr,/7T)"* deduced by this method is 
at variance with experimental results by a factor of two. (Here U denotes 
the speed at which the bubble rises, p the density of the liquid, 7, the 
equivalent radius (3 x volume/47)"® of the bubble, and 7 the surface 
tension.) 

Rayleigh solved the problem of the oscillation of a liquid or gaseous 
sphere at rest in an infinite medium otherwise at rest (see Lamb 1932, 
p. 473). ‘The more general solution for the oscillations of a gaseous sphere 
moving with a uniform velocity U in an infinite medium, to be obtained 
in the ensuing sections, reduces to Rayleigh’s result when U is equated 
to zero. 
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There has been a considerable amount of experimental investigation on 
the rates of rise of bubbles in various liquids (Hoefer 1913; Miyagi 1925; 
Bryn 1949; Datta et al. 1950; etc.). Rosenberg (1950) and particularly 
Haberman & Morton (1953) have conducted by far the most systematic 
experiments covering the spectrum of most of the pertinent variables. ‘The 
latter have also made an exhaustive study of the literature, collating the data 
of the previous experimenters and integrating these data with their own. 
Rosenberg improved the experimental technique and repeated the experi- 
ments on the rate of rise of air bubbles in water over a wide range of bubble 
sizes. He suggested the use of the three dimensionless parameters: drag 
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Figure 1. Drag coefficient as a function of Reynolds number for ellipsoidal and 
spherical-cap bubbles in various liquids. (Courtesy of David Taylor Model 
Basin, Navy Department.) 


coefficient C,, = (8/3)(gr,/U?), Reynolds number Re = 2pUr,/u, and a third 
parameter M = gu'/pT* for the description of bubble motion in liquids. 
Here g denotes the acceleration due to gravity and pw the coefficient of 
viscosity. Haberman & Morton, using the same experimental technique as 
Rosenberg, investigated the effect of variation of liquid properties on the 
motion of air bubbles and also made an experimental evaluation of the wall 
effect. They used eleven fluids with different properties (the same fluid at 
different temperatures being considered to be a new fluid) and three tanks of 
varying dimensions. In an attempt to correlate their data, they plotted their 
results in terms of the basic dimensionless parameters cited above. Figure 1 
gives the experimental results in terms of the drag coefficient, Reynolds 
number and M, whereas in figure 2, the data are presented in terms of the 
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drag coefficient, the Weber number and M. Examination of these curves 
indicates that there is no completely systematic arrangement with changes in 
M. Accordingly, it can be concluded that the variables considered in these 
non-dimensional parameters are insufficient for a complete description of 
bubble motion. Nevertheless, there is a definite trend towards categorizing 
the bubbles in groups according to the values of M. Figure 2 shows good 
correlation among the ‘fast’ fluids (MV = 10-! to 10-%) and sets a clear 
distinction between these fluids and the ‘slow’ fluids (WM = 10-7 to 10°). 
This paper is mainly concerned with the ‘fast’ fluids. 
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Figure 2. Drag coefficient as a function of Weber number for air bubbles rising at 


their terminal velocity in various liquids. (Courtesy of David Taylor Model 
Basin, Navy Department.) 


In examining the two plots, it is seen that the minimum drag attained 
by bubbles in various test liquids occurs at different Reynolds numbers, 
but at essentially the same Weber number (1-0 to 1-3). The correlation 
of all the curves for the ‘fast’ fluids above a Weber number of about 1-3 
indicates the importance of hydrodynamic and surface-tension forces in 
this range. 


2. EXPERIMENTS 


The experimental study was conducted to ascertain a conjecture of 
Professor von Karman that the Weber number should be a significant 
parameter involved in the stability of gas bubbles rising in liquids. 
Accordingly, measurements of the size and terminal velocities of bubbles 
which just become unstable were made in twelve liquids. In addition, the 
physical properties of the liquids were measured. 
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Figure 3 (plate 1) is a photograph of the experimental arrangement. 
The cylindrical tube is 5-7 cm in diameter, which is sufhcient to make 
any wall effect negligible. ‘The bubble-generating device consists of a 
7-section made of capillary tubing of uniform bore (0-7 mm inside diameter) 
and two rubber syringes. ‘This is a modification of the system employed by 
Saffman (1956). A shop microscope was used to measure the length (from 
which the equivalent radius is calculated) of small bubbles formed in the 
vertical section of the 7-shaped capillary. Fine black thread tied about 
the tube at two points a known distance apart (50-2 cm) provided the 
distance markings for the velocity measurements. ‘These were placed 
within the region where the bubbles had attained their terminal velocities. 
The time taken for the bubbles to traverse this distance was measured by 
means of a stop watch. ‘This measurement is the least accurate in the 
experiment. An estimate of the accuracy attainable is approximately 5°, 
for the more rapidly moving bubbles. ‘The viscosity of the liquids was 
measured by means of an Ostwald viscometer, and the surface tension by 
a DuNuoy tensiometer. ‘The density of each liquid at the appropriate 
temperature was obtained from the Handbook of Chemistry and Physics. 

The experimental procedure consisted first of measuring the physical 
properties of the test liquid. Air was bubbled freely through the liquid 
in order to reduce any tendency towards diffusion through the bubble 
surface. A single bubble was then formed in the vertical section of the 
capillary 7 by applying pressure to the syringe containing air until the 
desired quantity of air had penetrated the liquid in the vertical section 
of the 7, at which time the syringe containing the test liquid was squeezed, 
thus cutting off the desired volume of air. By continued pressure on the 
latter syringe, the bubble could be made to advance up the capillary tube 
to a point where its length was measured by the shop microscope (or a rule 
with fine divisions, for larger bubbles). Following this operation, the 
bubble was forced into the fluid medium at the orifice and carefully released 
from the orifice. The time taken for it to cross the markings was measured. 
Allowing sufficient time for the wake effects of the previous bubble to be 
dissipated, another bubble was formed and the process repeated. Starting 
with small bubbles and forming successively larger ones, a ‘critical’ size 
was found for each liquid at which the bubble would just begin to oscillate. 
Many bubbles formed about this size then made possible a better determina- 
tion of the critical conditions for each liquid. Upon completion of the test 
for a given liquid, a second temperature reading was taken, and, if any 
difference from the initial reading was detected, the physical properties of 
the liquids were obtained by interpolation from tables using the average 
temperature. 

From estimates of the experimental accuracy of the individual measure- 
ments, it is expected that probable errors of the order of 10% may occur in 
the calculation of the critical Reynolds and Weber numbers from the data. 
It is clear from this fact that only significant trends were being sought 
from these experiments. 
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Since it has been found that slight concentrations of impurities in 
liquids have profound effects on the motion of gas bubbles through them 
(Gorodetskaya 1949; Stuke 1952), great care was exercised in ascertaining 
the purity of the test liquids. Figure 4 compares the drag coefficients of 
bubbles in water containing slight concentrations of surface-active substances 
with those in filtered water. ‘Tap water and test liquids obtained from 
ordinary storage drums were found to behave in a manner analogous to 
water containing surface-active substances. ‘The drag coefficient in the 
impure liquids is seen to follow the curve for rigid spheres up to a Reynolds 
number of 200 or so. Of particular importance is the fact that the drag 
coefficient of bubbles in pure liquids is considerably smaller than that in 


the impure liquids. 
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Figure 4. Drag coefficient as a function of Reynolds number for bubbles rising at 
their terminal velocity in water containing various surface-active materials. 
(Courtesy of David Taylor Model Basin, Navy Department.) 


Aside from the criterion of purity, an attempt was made to select different 
liquids with some similar physical properties. For example, methyl alcohol 
and ethyl alcohol have essentially the same surface tension and density, 
but radically different viscosities. 

Since the experiments were designed to reveal the significance of the 
Weber number in determining the instability of bubbles, liquids of 
relatively low viscosity were selected to reduce the role played by the 
Reynolds number. Indeed, as discussed in a previous section, bubbles 
rising in very viscous liquids (roughly, MM > 10 *) do not become unstable. 

From the experimental determination of the critical bubble size and 
terminal velocity, as well as the properties of the various liquids, the critical 
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Reynolds, Weber, Bond r,(pg/T)'? and Froude U?/gr, numbers were 
calculated (table 1). The values obtained by other experimenters are also 
listed. Figure 5 is a plot of these results in terms of the Weber number 
and Reynolds number. ‘lhe average curve drawn through the experimental 
points is seen to give the strong indication that there are two distinct criteria 
for determining the onset of instability in gas bubbles rising in liquids of 
relatively low viscosity. ‘There is one branch of the curve that is independent 
of the Weber number, with the average value of 202 for the Reynolds number, 
whereas the other branch (at Re > 202) is independent of the Reynolds 
number and has the average value of 1:26 for the Weber number. 
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Figure 5. Stability curve (critical Weber number as a function of Reynolds number 
for various liquids). 


It seems significant that the former branch is made up not only of the 
more viscous liquids but also the impure ones. This may be taken as 
added evidence that impurities affect particularly the conditions at the surface 
of the bubble, introducing, it would seem, an increased apparent viscosity. 

When the nature of the stability curve became known, an experimental 
technique was developed to shed further light on the distinction between 
the two types of instabilities. Vegetable dye was added to the liquid in 
the syringe attached to the vertical section of the bubble-generating device 
(figure 3). By carefully forming a bubble, it was possible to see a closed 
wake behind a stable bubble as it rose. By closer examination, the dye could 
be seen to be rolled up into a vortex ring symmetrically situated with respect 
to the bubble. It is significant that the wake could be seen only in the impure 
and more viscous liquids. When the bubble size was increased to give 
critical conditions, the bubble oscillated and the wake dispersed. When 
this same experiment was attempted with distilled water, the following 
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sequence of events was seen to occur: (a) as the bubble left the orifice, 
the dye formed a closed wake; (4) after rising an inch or so, the bubble 
oscillated several times and shed the wake; (c) upon shedding the dye, it 
rose uniformly and rectilinearly; and (d) if the bubble size at the orifice 
was just below the critical size for distilled water, it would begin to oscillate 
before it had reached the surface. A possible explanation of this behaviour 
is as follows: (a) the presence of the dye in the otherwise pure liquid tended 
to transform its characteristics to those of the impure liquids; (4) in 
accelerating from rest at the orifice to the terminal velocity for distilled 
water, the bubble had to pass through the ‘ critical Reynolds number’ (202), 
and since it possessed the properties of an impure liquid, the bubble oscillated 
enough to shed its wake; (c) with the dye shaken from the surface of the 
bubble, the impurities were in effect removed, so that the bubble returned 
to its normal course in the pure water, i.e. it rose rectilinearly; and finally, 
(d) since the Weber number of a bubble increases (slightly) as it rises, owing 
to the reduction of the hydrostatic pressure and subsequent increase in 
volume (and therefore velocity), a bubble which was just sub-critical at the 
orifice might become critical as it rose. Fairly conclusive evidence regarding 
the action of the dye as an impurity is provided by the fact that if no dye 
was added to the liquid the short oscillation an inch or so above the orifice 
was never observed. ‘he experiment described above seems to support 
the existence of two different types of instability, both observable with one 
bubble. 

Referring back to figure 2, it is seen that the minimum drag for bubbles 
in liquids with small values of M occurs at a fairly universal value of the 
Weber number. Inaddition, the value of this Weber number is in agreement, 
within the limits of accuracy, with the experimentally determined critical 
Weber number for oscillations, i.e. 1-26. ‘The fact that the minimum drag 
coefficient and the onset of oscillation of the bubbles occur at almost the 
same Weber number is construed by the authors as providing additional 
evidence for the existence of a critical Weber number, since the oscillatory 
trajectory would contribute substantially to the observed rapid increase 
in drag with increasing Weber number beyond the critical value. Of course, 
the deformation of the bubbles at these Weber numbers ts also a contributing 
factor to the increase in drag. 

The stability curve (figure 5) is not quite analogous to the familiar 
stability diagrams in which a curve separates regions of stability from 
regions of instability. The actual meaning of the present stability curve 
is that for any liquid capable of producing unstable bubbles, the critical 
Weber number and critical Reynolds number must be such as to fit 
somewhere on this curve. Indeed, bubbles rising at Reynolds numbers less 
than 202 in Reynolds-number-dependent liquids are invariably stable, 
independent of the value of the Weber number, while bubbles 
rising in pure, relatively inviscid liquids at Reynolds numbers greater 
than 202 remain stable until they attain the value of the critical Weber 
number, 
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From table 1, it may be seen that there is no systematic variation of the 
Bond number with the Reynolds number at critical conditions, but that the 
critical Froude number increases fairly regularly with increasing Reynolds 
number. ‘The latter result may have been expected, since the Froude 
number is inversely proportional to the drag coefficient. ‘That the high 
values of the Foude number occur for liquids which have the Weber number 
as a stability criterion corroborates the contention that instability in these 
liquids begins at low drag coefficients (actually, almost minimum drag). 

The duality of cause for instability is not peculiar to bubbles rising in 
liquids. Similar results have been observed in investigations on the stability 
of liquid jets issuing from circular orifices (Richardson 1950, p. 100). 
It has been noted that falling drops break up at a critical value of the Weber 
number (Lane & Green 1956, pp. 186-192). 

The following conclusions are drawn from the present experimental 
results. 


1. A universal stability curve which determines the critical conditions 
necessary for the onset of oscillations of bubbles in most common liquids 
appears to exist. 


2. The character of this curve strongly suggests that there are two 
distinct criteria for instability : 
(a) a critical Reynolds number (202) for the impure and somewhat 
more viscous liquids ; 
(6) a critical Weber number (1:26) for pure, relatively inviscid 
liquids. 


3. ‘There is no exact method for determining which of the two criteria 
will prevail in a given liquid. For the pure liquids, a fairly good require- 
ment that the Weber number be dominant is that the M number should 
have a value less than 10-°. Unfortunately, the significant properties 
of ‘impure’ liquids do not seem to be reflected in their values of M, so 
that the liquids falling into this category can be identified only by 
experience. ‘lhis implies that the important properties affected by 
impurities are not simply viscosity and surface tension, at least as these are 
usually measured. 


3. "THEORY 


In this section the stability of a deformable sphere in an inviscid, 
incompressible, irrotational, uniform flow to a general small perturbation 
in shape is calculated. In applying the results to the bubble problem, it 
is assumed that the hydrodynamic pressure and surface tension represent 
the most significant effects. ‘lhe results indicate the possibility of small, 
high-frequency oscillations of the bubble at the lower Weber numbers, 
and divergence of the solution at Weber numbers above a critical value, 
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It will be argued that the region of small oscillations, when the damping 
effects of viscosity are taken into consideration, represents a stable motion 
of the bubble, whereas the divergent solutions actually imply the observed 
unstable motion. 

However, a subsequent investigation of the time-independent equations 
shows the existence of large deformations of the bubble prior to the onset 
of the unstable motion, which is not compatible with the approximation of 
perturbing an essentially spherical bubble. ‘Vhis deformation and_ its 
possible effects on the stability calculation are then determined by approxi- 
mate methods. From this, it is concluded that the deformation of the bubbles 
serves to introduce quantitative, but not qualitative, changes in the stability 
calculation. 

This approach to the problem differs significantly from that employed 
by Saffman (1956) in two respects. Firstly, it is assumed that the flow about 
the entire sphere is essentially inviscid, and secondly, in addition to that 
mode which represents a lateral velocity perturbation, changes in shape of 
the bubble surface are considered to be important. 


A. Theoretical assumptions and development of the boundary condition 

The theoretical analysis is concerned only with the instability of bubbles 
due to hydrodynamic pressure and surface-tension effects. ‘The correlation 
obtained for the critical conditions in the pure liquids with less than 10~% 
in terms of the Weber number (figure 5) gives significant evidence of the 
primary roles played by the hydrodynamic pressure and surface tension 
in the instability. ‘That viscous effects are relatively unimportant in this 
phenomenon is, at first sight, somewhat surprising, considering the low 
Reynolds numbers encountered. However, the drag coefhcient of the 
bubbles just prior to the onset of oscillations is approximately one-fifth 
that of a rigid sphere at the same Reynolds number. In fact, rigid spheres 
do not attain such low drag coefficients until Reynolds numbers of the 
order of 4x 10° are reached. At these Reynolcs numbers the pressure 
distribution about the sphere is essentially that obtained from potential 
theory. On the basis of this evidence, the theeretical model will have the 
form of a deformable sphere moving uniformly in an inviscid, incompressible 
fluid of infinite extent. 

Gravitational forces and the pressure of the gas within the sphere due 
to perturbations in the shape of the sphere will be neglected. ‘Their effects 
may be shown to be much smaller than those of surface tension and 
hydrodynamic pressure for the size bubbles considered (r, < 0:09 cm). 

Since the bubbles do not change volume appreciably as they rise, the 
perturbation of shape will be performed so as to maintain constant volume. 

Employing the usual spherical coordinates with the origin situated at 
the centre of the bubble, the equation for the surface of the bubble is denoted 


. 


by 


R(0, 4,1) = Ry+ > A,(t)S,(0, 4), (1) 


n=1 
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where Ry is the equilibrium radius of the bubble, 4,,(¢) are small time- 
dependent quantities such that their squares may be neglected, and the 
S,,(@,) are spherical surface harmonics represented by 


S,,(4,¢) = cy P,(cos#)+ > (c,, cos md + d,, sin md) P'(cos @). (2) 


In order to maintain constant volume A,(t) is taken to be zero; thus the 
sum in equation (1) starts with n = 1. Let 

A(t) = a,,e” (3) 
where the coefficients @,, are small (compared to R,) complex quantities, and 
A is a complex number to be determined. If the real part of A is negative, 
the ensuing motion is stable, whereas positive values of A will indicate a 
divergence of the mode in question. ‘lhe surface of the perturbed bubble 
may accordingly be represented by the equation 


F(R, 0,4, t)=Ry+e* > > a 


[m= 


nm o. cosmo — R = 0, (4) 
where it is understood that the argument of the P” is cos@ throughout 
the analysis. ‘lhe boundary condition at the surface of the bubble may be 
derived from the condition 


DF/Dt = 0, (5) 


where D/Dt is the usual convective derivative. Since we are assuming 
potential flow about the sphere, the velocity components at the surface of 
the perturbed sphere are of the forms 


v, = 3Usiné+v,+O(a.,,), v, = v. + O(a;,,,), Ug = Us t+ O(Gim), (4) 


where the primes denote quantities of the first order in the perturbation 
amplitudes a,,,. If equations (4) and (6) are substituted into (5), the 
boundary condition on the surface of the perturbed bubble is obtained, to 
first order, as 


ve) n 
[Vela = Ae > > Gam Pr cosmo + 


n=1 m=U 


3 oF } m-+- | n+1)\(n+m 
n 5 ny De eee 7 ) p", - oP m) pu _ Joos md, (7) 
where use has been made of the recurrence relation for sin@dP’"/dé (Morse 
& Feshbach 1953, p. 1326). 

In the following section, a potential will be designed to satisfy the 
boundary condition given by (7), and from this potential the hydrodynamic 
pressure will be calculated. 


B. Development of the solution 
The velocity potential must satisfy V°2 = 0. ‘The potential may be 
written as 


(p }U(2r+ R,/r?)cos6+ 0", (9) 
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where ®’ contains terms of the first order only, and the first term is the 
familiar velocity potential for the flow about a sphere of radius Ry. Aside 
from the first-order terms in (7), it is seen that the first term in (8) also 
contributes a first-order term in the radial component of the velocity when 
it is evaluated on the surface of the — sphere, for 


[v, lz = lv,)p,4 | | AR = -—3 a cos 6 e y y a, P“ cosmd. (9) 
Ro 0 = ’ 


‘herefore, in addition to finding a potential aes will = equation (7), 
the potential must also cancel the contribution indicated in (9). If use is 
made of the recurrence relation for cos 6P", equation (9) may be rewritten 
l ae (n—m+1) (n+m) 
an ae beta “ ‘ Gs UR Pp" L | em ie d. 
Vr-lr = —35 ¢ a, : , |cosmd. 
[ rr Ry 2 * nm (2n ni 1) (2n +1) ‘s 1) | 
(10) 
If one subtracts (10) from (7), an equivalent boundary condition remains 


to be satisfied, viz. 











ae 3 l —_—-! 
[o]p=Ae*> Sd a,,, P? cosmo+5—e" > DY an» 
n-lm=0 2 Ry n=1m=0 
(n+2)(n—m-+1) (x —1)(n+m) 
x Pp" |.- ——— P_, |cos md. 1] 
| 2n +1 | 2n+ 1 sili » a 
The appropriate velocity potential is seen to be 
(U R: . 2 é£ 
O= —-< 2r+ — }cos@+rAe* YY —“ —— P" cos md + 
2 ( ” :) nel = (a v Wea " 
$U iS 
ris: A? s 
( a 
2M scinco 
(n—m+1) Rv (n—1)(n+m) 
x} ——— SS Pr" 1- ”) Re tp" | |cosmd}. (12 
| n+) ett’ *t! “atte 8 8 A ata) 
The pressure is calculated from 
p/p +4v?+e@/ot = constant. (13) 


Inasmuch as a first-order theory is being constructed, it is clear that the 
v, will not contribute to the fluid pressure, 


velocity components 7, and 
since each is already of the first order in a,,,,._ ‘The velocity component v,, 


on the other hand, is, from (12), 


= =f sin 6 — 3 nay Ly s a,,, sin@P" cosmd ? = (14) 
2 21 o n= 0 R 6? 
so that 
* Lean =- : U? sin®6 - 
; e Ye > y a, sin?¢P” cos md - si sin Pil (15) 
4 Ry n=1 m=0 nee 2 . " aa : 


where terms of the first order only have been retained. If use is made of 
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the recurrence relation for sin?@P’ (Morse & Feshbach 1953, p. 1326), 
it can be shown, after some algebra, that the fluid pressure as calculated 
from (13) is 


p to Be, {[(nt+m-+ 1)(n+m-+ 2) 
Cl wits aw ees oS Ss pb es oat Bd A Ms ae 
| - p 2 8 U in fe r 4 a * Se Sam (2n + 1)(2n4 3) T 


p 
a (n—m—1)\(n—m) (n- L)(n + 2)(m +m +1) 


-(2n—1)(2n S | i Wa { 1)(2n 3) 


_ (n-t 1)*(n* se m*) p 
n(2n + 1)(2n—1) 


eo ee) a 
1)(2n + 3) “a + 1)(2n +3) m 
k -m)\(nt+m—1) (n—1)(n m)(n + mt 


ban 2] ge “Ae Id; 
(2n + 1)(2n— 1) (2n + 1)(2n—1) ‘ee 


| 


oe: (n+m) n(n —m-+- 1) 
Unde"! = > th E — P™ = ees P_, |cos me + 


+ 
> 
t 
& 
> 
= 
SS 
~ 
VU 
a 
~ 
7 
= 
= 
= 
= 
oO. 
= 


: oe PP 1- pe - Joos mp, (16) 


2n+1 sie n(2n + 1) 


where the terms 0®/ct and d@’/d(cos @) were calculated from (12). 
The action of surface tension at an interface between two fluids is to 
cause a discontinuity in the pressure across the surface, given by the relation, 


p:—p = T(1/R, +1/Ry), (17) 


where 7' is the surface tension, and R, and R, are the principal radii of 
curvature, considered positive if the centres of curvature lie on the side to 
which the subscript 7 refers. For the perturbed sphere in the present 
problem, the following form for the sum of the curvatures, correct to first 
order, has been derived in Lamb (1932, p. 474), 
| Z . 2 (n—1)(n+2) 
~+—=75 +e7> Y —— Gan Pi cos md. (18) 
R, R, Re, a repel Rk? 
Substituting (16) and (18) into equation (17), one has for the equilibrium 
of the pressure terms that are independent of time, 
oU2 9 27 
i os es Retntis 2 ciel sae ( 
P, =P. i) pl sin“@ Pi a0; in =O TA me Olen =O? (19) 
8 RG 
where the last two terms refer to the first-order fluid pressure and surtace- 
tension stresses which are time-independent and symmetric with respect 
to the direction of motion. Equation (19) permits calculation of the 
deformation of the bubble as a function of the Weber number, which will 
be considered in a later section. For the equilibrium of the time-dependent 
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pressures terms at the surface of the bubble, one has 


, pee (n-+m+2)\(n+m+1)\n 3 , (7 + m+ 1)(2n +1) 
oof Oyen ee = WY’ Gy 1m eee 
+ “> (2n + 3)(2n +5) 2 “(n+ 1)(2n + 3) 
Y We C(n+m+1)\(n+m+2) (n-m—1)(n—m) 
~ Wa, | —as ——.-—- on 
4 si | (2n + 1)(2n + 3) (2n—1)(2n+ 1) 
(n—m+1)(n+2)(n+m-+ 1) (n — 1)?(n? — m?) 
(2n + 1)(2n +3) n(2n —1)(2n + 1) 


De. | 
{ = 4 iMa+t? 
+ G,, ol 5 jz t(e-le a 


(1—m) 9 W2 (n—m—1)(n—m)\(n 
ens - = 
+ gece (2n — 1)(2n—3) 


4 


o 
5 WA ay—1 


0 (20) 
mh = Ps 

where the conditions of orthogonality of the associated Legendre functions 
and the trigonometric functions have been employed, A’ -= A(pR})!?/T'?, 
and W denotes the Weber number. ‘The ranges of n and m are n } and 
0 <m <n, and the a, for n < 0 are considered to be zero. ‘That this 
equation should be equated to zero follows from the neglect of the perturba- 
tion pressures of the gas within the bubble. for the case of zero forward 


velocity (U = 0), equation (20) reduces to the result obtained by Rayleigh, 
Viz. ( | 1\(n+2 
: (n+ 1)(n—1)(n+2),, 
A =. « = — Th een l ° (21) 
oR 


I 
It is clear that for this case there is no coupling of the modes, or otherwise 
stated, the coordinate system selected for the problem consists of the normal 
coordinates. ‘lherefore, it is seen that the effect of a uniform translation of 
the bubble is, effectively, to couple the possible modes of oscillation, or at 
least the modes dependent upon 6 (i.e., the indices 2). In seeking solutions 
to (20) for the stability parameter A’ as a function of the Weber number, 
it is clear that there will be distinct solutions A’ for each m, with.all the 
coupling modes (”) contributing to these solutions. ‘Therefore, the 
technique will be to assign a value to m and to solve numerically the resulting 
infinite determinant. ‘lhe analysis will be concerned mainly with the smaller 
values of mn and m, because these deformations most nearly simulate those 
observed in experiments, and also because the surface tension, which 
maintains stability of the shape of the bubble, is much more dominant than 
the hydrodynamic pressure for the higher-order deformations (m large), 
as evidenced by the dependence on rn? of the term due to surface tension in 
(20). Actually the modes x = 1, m = 1 and » = 2, m = 1 are most closely 
representative of the observed motion of the bubble once instability occurs. 


C. Numerical solution of equation (20) for asymmetric modes. 

It has not been possible to obtain a general solution of (20) in closed 
form. Indeed, there is no general theory available for the numerical 
solution of a‘five-term recursion relation. However, for the neutral point 
of the stability calculation A’ = 0, equation (20) reduces to a three-term 
recursion formula involving the Weber number and the a,,,,.. ‘Uhis relation 
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is analogous to the three-term recursion formula obtained in the calculation 
of the stability of the solutions of the Mathieu equation (Morse & Feshbach 
1953, p. 557). ‘There it is shown that the a,,,, converge only for certain 
values of W* and ay,/a,,. ‘These characteristic values are found by the 
method of continued fractions. ‘This technique is the same as forming 
the determinantal equation from (20) with A’ = 0 and solving the resulting 
polynomial in W? for the characteristic values. Of course, the determinant 
contains an infinite number of terms for each value of m, but it is found that 
rapid convergence for the characteristic values W? is obtained by cutting 
off the determinant in successively larger segments. ‘lhis process of cutting 
off the determinant is equivalent to terminating the continued fraction. 
‘he rapid convergence is to be attributed to the term (m—1)(m +2) in (20) 
which predominates over the other coefficients which vary as n for large 
values of for any given m. ‘Uhe results of this calculation give values of 
the critical Weber number for each m. ‘To see that values of the Weber 
number smaller than this lead to stable solutions, the general equation (20) 
is solved by an analogous method to obtain a curve of A? vs W?. It is 
found in expanding the determinants that only even powers of 4’ occur, and 
further, regardless of the mode selected or order of the determinant included 
in the calculation, the values of A" are always real. ‘That this should occur 
is suggested by the absence of any damping terms in the problem. From 
these facts it is clear that all solutions are either purely oscillatory (A’? < 0) 
or purely divergent (A? > 0). 

‘The critical Weber number calculated as described above, tor m = 1, 
is 1-65. Figure 6 shows how A’? varies with the Weber number as calculated 
from both the 2x 2 and 4x4 determinants for m= 1. Going to higher 
order determinants gives results only negligibly different from those obtained 
from the 4 x 4 determinent. ‘lhe amplitudes converge exceedingly rapidly. 

‘The curves of figure 6 indicate the possibility of stable oscillations 
of the bubble in the lower modes for Weber numbers ranging from zero to 
critical. Of course, these oscillations are not observed in the experiments, 
since in this range the bubbles are seen to rise rectilinearly. ‘This discre- 
pancy may be explained by the absence of viscous effects in the theory. 
For Weber numbers of the order of unity to critical, the predicted values of 
A” vary from —0-5 to —O-1, Since the non-dimensionalizing factor in 
(i.e. 7'/pR}) is of the order of 1-5 x 10°, the actual values of 7A range from 
250 to 120 approximately. ‘This results in periods of oscillation between 
0-025 sec and 0-05 sec, which corresponds to the time required for the 
average bubble to travel only four and eight bubble diameters, respectively. 
At these relatively high frequencies of oscillation, it may be expected that 
viscosity would damp the motion rather completely, since the theory deve- 
loped by Rayleigh (see Lamb 1932, p. 640) may sensibly be applied. 
Consequently, the purely imaginary values of \’ encountered in the analysis 
are to be considered as representing stable motion of the bubble, with 
the possible exception of the points just short of the critical Weber number 
where A’ becomes equal to about 0-V1lz. ‘Then, perhaps, viscosity would 
be less significant in damping the stable oscillations. However, since the 
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objective of the theoretical analysis is to calculate the Weber number at 
which the bubbles become unstable, it does not seem necessary to dwell 
on the question of whether the observed instability occurs for A’? equal to 

0-01 or +0-01. For Weber numbers above critical, the motion is a 
purely divergent, asymmetric one which presumably represents a darting 
motion of the bubble to one side of its original path. Motion pictures of 
unstable bubbles leaving an orifice with vegetable dye in their wakes indicate 
that the amplitude of the oscillation is several bubble diameters. ‘This 
is indicative of a non-linear motion, or at least, a motion that could not be 
derived from small-perturbation techniques. 
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Figure 6. ‘Uheoretical stability curve. 

‘lo summarize the preceding paragraph, the predicted high-frequency 
steady oscillations of the bubble at small Weber numbers is interpreted to 
correspond to the steady rectilinear motion actually observed. ‘The 
divergent motion predicted for Weber numbers above the critical value 
(or, what is equivalent, the slow steady oscillation predicted just below the 
critical W) is interpreted to represent the large-amplitude oscillatory motion 
actually observed in that regime. 

In closing this section of the analysis, it is necessary to mention once 
again that only’ very slight steady-state deformation of the bubble was 
considered, consistent with the neglecting of the squares of the perturbation 
amplitudes. In practice, it is found that the sub-critical bubbles are some- 
what more deformed in the shape of oblate spheroids than this calculation 
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may legitimately permit. ‘This question will be considered in more 
detail in the following section. However, it may be said that the analysis 
presented here serves to explain qualitatively a possible mechanism for the 
instability of bubbles under the basic assumptions considered. 


4. DEFORMATION OF BUBBLES AND ITS EFFECTS ON THE STABILITY 

The deformation of bubbles from the spherical shape as they rise 
uniformly in an inviscid fluid under the action of the hydrodynamic pressure 
and surface tension may be calculated from the steady-state terms of the 
previous calculation (see (19)), at least to first order in the perturbation 


O6rf 











Figure 7. Bubble deformation calculated by approximate methods. 


amplitudes. ‘The relations for pi; — .. ,, —9) and 7; — 9, » —o) may be calculated 
from (16). If the resulting system of linear, inhomogeneous algebraic 
equations is solved by Cramer’s rule, one obtains 
ay . 3 W(3-58W?—-18) 22 
R,  4(5:52W4—43-6W24+72) (32) 
approximately. Again the infinite set of algebraic equations has been 
terminated and rapid convergence has been obtained. ‘The ratio dg9/Ry 
is a measure of the deformation of the sphere into an ellipsoid at various 
Weber numbers. Equation (22) is plotted in figure 7 (dashed line). ‘That 
such large deformations in shape are predicted by the linear theory is, at 
first thought, quite alarming, for it places in question the application of the 
stability calculation of the previous section to gas bubbles rising in liquids. 
It is observed that just sub-critical bubbles are quite deformed in the 
experiments—almost as much as indicated by (22). Therefore the possible 
effect of ellipticity of the bubbles on the quantitative values of the stability 
criterion of the asymmetric modes will be dealt with in this section. It will 
be shown that the qualitative behaviour is not altered. 
To obtain an estimate of higher-ordered effects in the deformation of 
a bubble in an inviscid fluid, one may start with the pressure distribution 
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about an oblate spheroid in such a fluid, and knowing the principal radii 
of curvature at the stagnation point (6 = 0) and equator (6 = 7/2), simply 
require that the general equilibrium condition (17) be satisfied only at these 
points. ‘lhe hydrodynamic pressure at the equator is given by 


p = \pU2—IpU%1 +h.) (23) 
where k, = (0-622E — 0-122) (23 a) 
and E is the ratio of semi-major to semi-minor axes (b/a). ‘The values of k, 
are tabulated in a paper by Zahm (1926). ‘The sum of the reciprocals of 
the principal radii of curvature at the equator is (6/a?+1/b), while that 
at the stagnation point is 2a/b*._ Therefore, for equilibrium, 

p; = 3p U? + 2aT |b? (24) 
and p; = pU?—3pU2(1+k,)? + T(b/a? + 1/5), (24a) 
and since the gas pressure within the bubble (;) is uniform, 

T(2a/b? — b/a? —1/b) = —3pU%1+k,)?. (25) 
‘To maintain constant volume, we have 
al? = Bi, i.e. b = R, E"*, (26) 
The Weber number obtained by solving (25) and (26) is 
2 (E?+1-2/E) 


~ Baek” 


Ww? 


which gives the deformation of the bubbles as a function of the Weber 
number. Since the semi-minor axis of the spheroid is equal to R,/E** 
and the amplitude ratio ay /Ry =(a/R)—1), we may use the relation 
Ay)/Ry = (E-*3—1) in conjunction with (27) to plot this latter equation 
adjacent to (22) in figure 7 (solid curve). ‘The results show that (27) verifies 
the linearized theory for small Weber numbers, but gives smaller 
deformations for somewhat larger values. 

The existence of distortion in the shape of the bubbles would make the 
perturbation of an oblate spheroid, rather than a sphere, seem a more 
reasonable attack on the problem of the stability of the bubbles. ‘The 
equations, employing oblate spheroidal coordinates, have been derived. 
Though not amenable to complete solution immediately, this analysis further 
verifies the qualitative correctness of the theory developed for the sphere, 
and together with the ‘two-point’ calculation mentioned above, suggests 
a method for deriving better quantitative results from this theory. Briefly, 
this method consists in replacing the form (3/2)Usin@ for the tangential 
velocity by the increased value appropriate to an oblate spheroid of any given 
eccentricity (Zahm 1926), in the boundary condition (5). If this is done, 
it is found that instead of having 9/4W7? equal to the characteristic value of 
(20) for the neutral stability point (A’ = 0), the square of the increased value 
replaces the coefficient 9/4. ‘Thus, an improved estimate of critical Weber 
number may be obtained for any assumed deformation E (equal to the ratio 
of semi-major to semi-minor axes of the spheroid). But the appropriate 
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deformation is actually a function of Weber number, and may be estimated 
by the ‘two-point’ calculation. ‘lhe consistent case is the one for which 
the critical Weber number is just that required to produce the assumed 
deformation. ‘This turns out to be F = 2-20, for which W,, = 1:23. ‘This 
Weber number is in excellent agreement with the experimental value. 
Moreover, the deformation of bubbles in most of the pure, fast liquids just 
prior to instability, as estimated from photographs supplied by Haberman 
& Morton, is approximately given by & = 2:1. An outstanding exception 
is provided by bubbles in the fi/tered water employed by these experimenters, 
where F= 1-5. ‘The present authors employing double-distilled water 
available commercially, obtained a deformation E = 1-75. 

It is suggested from the results of these approximate calculations that 
the deformation of the bubbles prior to oscillation reduces the value of the 
critical Weber number as calculated from the perturbation of an essentially 
spherical bubble, thus bringing it into better agreement with observations. 
Throughout this section it has been implied that the Weber-number effects 
are largely responsible for the actual deformation of the bubble; hence the 
present agreement with experimental results supports all previous experi- 
mental evidence in establishing the dominance of these effects for bubbles 
rising in the pure, fast liquids at Weber numbers near and somewhat 
beyond the critical value. 


5. CONCLUSIONS 

The conclusions derived from the experiments were listed at the end of 
$2. From the theoretical analysis, the following conclusions may be drawn. 

1. The results derived from the investigation of the stability of a 
deformable sphere in a uniform, inviscid, incompressible and irrotational 
flow by the perturbation of the shape of the sphere, indicate, at least 
qualitatively, that the interaction of surface tension and hydrodynamic 
pressure provides a possible mechanism of instability. In fact, in any 
situation where the deformable sphere is not subjected to considerable 
changes in shape just prior to the onset of instability, this analysis would 
be — to give accurate quantitative results as well. 

‘The agreement between the approximate theory and the experimental 
ies for the steady-state deformation of the bubbles under the action of 
only hydrodynamic pressure and surface tension, for moderate values of 
the Weber number, provides additional evidence of the important roles 
— by these effects in the motion of bubbles in that range. 

The inclusion of the effects of deformation of bubbles by approximate 
cna in the stability calculation of the nearly spherical bubble yields 
a critical Weber number of 1-23, which agrees well with the experimental 


value. 
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NOTE ADDED IN PROOF 

In discussion of the stability curve (figure 5) since submission of this 
paper, the question has been raised whether the horizontal (IV = 1-26) 
and vertical (Re = 202) branches of this curve should not be extended 
to the left and upward, respectively, from their point of intersection. 
In answer to this question, the authors point out that neither they nor 
other experimenters have observed oscillations in these regions of the 
diagram. Asa bubble rises from rest in a liquid, its instantaneous Weber 
and Reynolds numbers trace a curve that rises toward the right from the 
origin in figure 5, its steepness depending on the fluid properties. For 
the experimental points plotted in figure 5, a Reynolds-number-dependent 
fluid produces instability when this curve crosses Re = 202; a Weber- 
number-dependent fluid does not cause instability as the curve crosses 
this branch, but only where it reaches WW’ = 1:26. ‘The ‘ more viscous 
liquids’ mentioned on page 28 seem to be those whose trajectory curves 
rise to the left of the point (202, 1-26) in figure 5, and, as mentioned, these 
do not exhibit instability at all. ‘To be sure, there are only a few data on 
these liquids; nevertheless, the extrapolation of the horizontal and vertical 
branches of the stability curve must be considered unwarranted at present. 
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SUMMARY 

When plane pressure waves in a duct reach an open end, they 
establish a complicated three-dimensional wave pattern in the 
vicinity of the exit which tends to readjust the exit pressure to its 
steady-flow level. ‘This adjustment process is continually modi- 
fied by further incident waves, so that the effective instantaneous 
boundary conditions which determine the reflected wave depend 
on the flow history. In the analysis of a nonsteady-flow problem 
by means of a wave diagram, it has been customary to assume that 
the steady-flow boundary conditions are instantaneously establi- 
shed. While this simplifying assumption appears reasonable, the 
resulting errors have been undetermined. It is the purpose of 
the present investigation to obtain improved boundary conditions. 
The results of a previous study of the reflection of shock waves 
from an open end have now been extended to other ‘waves of 
finite amplitude. ‘The reflected waves computed by means of the 
new procedure are in good agreement with experimental data 
observed in a shock tube for a variety of flow conditions. ‘The 
pressure variations in a reflected wave lag behind those derived 
in the conventional manner by the time in which a sound wave 
travels about one or two duct diameters. Such lags are small, 
but may occasionally become significant. As a consequence of 
the lag, certain discontinuities of the incident wave do not reappear 
in the reflected wave. ‘This improved understanding of the 
reflection process has made it possible to clarify some previously 
unexplained experimental observations. 


1. INTRODUCTION 


The reflection of pressure waves of finite amplitude from an open end 
of a duct must frequently be determined when quasi-one-dimensional 
nonsteady flows are being analysed by the method of characteristics. Such 
reflections are customarily computed under the assumption that the 
boundary conditions in nonsteady flow are the same as in steady flow 
(e.g. Shapiro 1954, p.960; Rudinger 1955a, p.59). ‘This assumption 
leads to convenient computing procedures, but represents only an approxi- 
mation to the actual phenomena. ihe incident waves change the 
steady-flow exit pressure which, through the mechanism of a complicated 
three-dimensional wave pattern, can then re-establish itself only gradually at 
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a rate that depends on the time for waves to travel across the end section of 
the duct. Since this adjustment process is continually modified by the 
incident waves, the instantaneous boundary conditions depend on the 
history of the flow. Application of the steady-flow boundary conditions 
implies that the adjustment is instantaneous, so that the deviations of the 
exit pressure from its steady-flow level are neglected. ‘lhe present investi- 
gation was undertaken to obtain improved boundary conditions which may 
be used instead of the conventional procedure or, at least, may permit an 
evaluation of the lag in the establishment of the steady-flow boundary 
conditions. Since wave-diagram procedures are based on the model of a 
one-dimensional flow, it is only necessary to find ‘effective’ boundary 
conditions. Details of the highly complicated three-dimensional flow 
phenomena in the immediate vicinity of the end section of the duct are not 
needed. 


‘The work started with an analysis of the refl 


lection of shock waves from 
an open end (Rudinger 1955 b), based on acoustic theory, and the results 
were in good agreement with experimental observations in spite of the 
finite amplitude of the waves considered. ‘lherefore, it seemed promising 
to extend the analysis to arbitrary incident waves. ‘The ‘ effective’ pressure 
variations that occur in the end section of the duct during the adjustment 
process following the arrival of an incident shock wave may be considered 
as the ‘ response to a step function’, from which the response to an arbitrary 
wave can be computed by means of Duhamel’s integral. It is implied in 
this analysis that the steady-flow boundary conditions are represented by 
a constant pressure at the duct exit equal to that of the surrounding gas, but 
this boundary condition applies only to subsonic outflow from the duct. 
(The restriction to subsonic flow must be made because no reflected wave 
can propagate into the duct if the outflow velocity is sonic or supersonic.) 
If the direction of the velocity is reversed, and the gas flows from the 
surrounding gas into the duct, the pressure at the end is not constant, but 
becomes a function of the flow velocity, even in steady flow. ‘The formation 
of a vena contracta, caused by flow separation at the inlet edge of the duct, 
further complicates the situation. ‘The computing procedure derived 
from the results of the acoustic theory with the aid of Duhamel’s integral 
strictly applies, therefore, only to outflow; but certain empirical modi- 
fications make it possible to deal also with inflow. 

An experimental check on the results is obtained by measuring the 
pressure variations of the incident and reflected waves at some distance 
from the end of the duct where the flow can be considered as one-dimen- 
sional. From the data for the incident wave alone, the reflected wave can 
be computed on the basis of the new and the conventional procedures, and 
the results compared with the experimental observations. 

2. DERIVATION OF THE COMPUTING PROCEDURE FOR SUBSONIC OUTFLOW 

In the previous study of the reflection of a shock wave from an open 
end (Rudinger 1955b), the incident wave was expressed in terms of its 
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frequency spectrum by the Fourier integral of a step function. ‘The 
Fourier integral of the reflected wave was then obtained from the known 
frequency dependence of the acousticimpedance of anopenend. Numerical 
evaluation of this integral yielded the effective variations of the exit pressure, 
and the results were in good agreement with experimental observations. 
For weak shock waves, there are noticeable differences between the observed 
reflected wave and that computed by means of the conventional procedure. 
These differences become insignificant for stronger shock waves (i.e. with 
pressure ratios larger than about two in air). This is a fortunate finding, 
because the use of acoustic theory implies isentropic changes of state and 
applies, therefore, only to weak shock waves. 

The deviations of the instantaneous pressure at the end of the duct 
p.(z) from the steady-flow value py were expressed in terms of a dimension- 
less time 

tT = at/D, (1) 
where ¢ is time measured from the instant of arrival of the shock wave at 
the exit, and a, is the initial speed of sound; the duct diameter D is the 
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Figure 1. Plot of the functions J(7+) and P(r). 


significant dimension that controls the rate at which the steady-flow bound- 
ary conditions are approached (see §3 for further remarks on this). The 
result was obtained in the form 


PAT)—Po = (Pi-Po)A(7), (2) 
where p;—p, is the pressure rise across the incident shock wave, and the 
time dependence J(7) is shown in figure 1 (see also table 1). Equation (2) 
represents the instantaneous boundary conditions which can be applied 
in a wave diagram. 

One may consider the function J(7) as the indicial admittance of the 
system, that is, the response of the exit pressure to an incident wave in the 
form of a unit-pressure step. The response to an arbitrary incident wave, 
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described at the exit by p(r)—p,) = F(r), can then be derived by means of 
Duhamel’s integral which may be expressed in the form (see, for instance, 
Karman & Biot 1940, p. 403) 


‘ “tT dF(6 A 
PAT) Dy F(t 9)I(7 — 7) + | ; : s ) (76) dé, (3) 





where 7, is the instant at which the incident wave begins to arrive at the 
end of the duct, and @ denotes the integration variable. ‘The wave diagram 
is started from an initial state of rest or steady flow; thus, the condition 
F(r) = 0 for + <7, eliminates the first term on the right side of equation (3). 
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Figure 2. Wave reflections from an open end. I and II indicate the location of 
pressure transducers. 


For 7>7,, the function F(r) cannot, in general, be prescribed in 
advance, but becomes available only as the construction of the wave diagram 
proceeds. Figure 2 shows a sketch of the wave diagram which includes 
several characteristics of the’ incident wave (P), P,, P...) and their 
reflections from the open end (Qo, Q,, Q,...), which propagate with the 
velocities u,, +a, and u, —a, respectively (x = 0, 1, 2...). The Riemann 
variables P,, and Q,,, defined by 

P,, = 2a,/(y—1)+u,, (4a) 
O, = 2a,/|(y— 1)—u,, (4 b) 
where y has its usual meaning, remain constant along their respective 
characteristics if the flow is isentropic and the duct has a constant cross 
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section. (‘lhese assumptions do not restrict the present analysis, which is 
required to find only the wave that reaches the duct exit regardless of the 
manner in which the wave is generated. See. $3, where nonisentropic 
flow into the duct is considered.) 

Regions between any two consecutive characteristics of the incident 
wave may be considered as small wave elements for which the derivative 
in the integrand of equation (3) can be replaced by the finite-ditference 
ratio AF/Ar. ‘The flow conditions remain constant along a characteristic 
of the incident wave only until the latter interacts with the reflected wave. 
In the interaction region, the strength of a wave element must be measured 
along a crossing characteristic (see Rudinger 1955a, pp. 31 & 32). ‘The 
strength remains constant if expressed as a change of the speed of sound, 
but not if expressed as a change of pressure. Instead of calculating the 
pressure changes in the incident wave as modified by the reflected wave, 
it is therefore preferable to derive the boundary conditions directly in 
terms of the speed of sound, so that /(7) must describe the incident wave 
in terms of this variable. ‘The strength of an element of the incident wave 
Aa, follows, from the property of characteristics (equations (4) with 
O = const.), as 

Aa, = Hy—1)(P, - P,-1)- (5) 
This definition, apparently, neglects the last piece of the characteristics 
before they reach the exit. In figure 2, for example, 

Ady = Ag9— A, 9 = 421-4, 
and the change of a between points 2, 1 and 2 is not considered. Actually, 
the incident wave element should be taken as a,—a,, where the prime 
indicates the flow conditions that would be established if the duct did not 
terminate at this location after the arrival of the characteristic P,. ‘There 


would then be no reflected wave between 1 and 2, so that O5 = O, = Op,, 
and, therefore a, = a),. This reasoning justifies the foregoing definition 
of Aa,,. 


Interaction of the incident and reflected waves also modifies the times 
at which the characteristics of the former arrive at the open end. ‘The 
actual arrival times, and therefore the values of Av, can be taken directly 
from the wave diagram (see figure 2). 

‘lo derive a computing procedure in terms of the speed of sound rather 
than pressure requires that the indicial admittance be expressed in terms 
of this variable. Equation (2) defines /(7) as a pressure ratio which is not 
exactly equal to the corresponding ratio of the sound speeds. Since the 
changes of state may be considered as isentropic, equation (2) may be 
rewritten in the form 

-pde)= Po _ (4 a)’"Y—-1  ar)- ao 


= Pi- Po 7 (a; Ay)?” ioe sl 








p, (6) 


which defines a correction factor 7s that is always less then unity, because 
ad; <a,7) <a. Immediately after the arrival of the shock wave, when 


a a;, one obtains 4 =1. The minimum value of % is reached when 
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a, approaches ay; but this error is insignificant because J(7) then approaches 
zero. It seems reasonable, therefore, that any significant error would 





appear near the middle of the range where a, = }(a;+4)). For this value, 
a series expansion of ¢# yields 
1/y+1\a;— a, ~ 
sh = 1- =| —— Bina / 
? Hy i) a. (7) 


For increments of a;— a, not larger than about 0-02a , the value of ¢ is then 
greater than 0-9 for any y larger than 1-1. Such errors of the indicial 
admittance curve (see figure 1) could not be detected experimentally. 
Throughout the following, the relation 


aT) — ay 














I(r) = a;— Ay (5) 
will therefore be used. 

— 
| T I(r) (7) = I(r) P(r) 
| ih ee ool 
| 0 1-000 0 1-6 (0-126 0-774 | 
| 0-1 0-938 (0-097 1:8 0-094 0-795 | 
| 0-2 0-870 0-187 2-0 0-070 0-811 | 
| 0:3 0-798 0-271 2:2 0-052 0-823 | 
| O-4 0-721 0-346 2-4 0-039 0-832 | 

0-5 0-642 0-414 2-6 0-029 0-839 | 
| 0-6 (0-570 0-475 2:8 0-021 0-844 | 
| 0-8 0-431 0-574 3-0 0-014 0-847 | 
EO) 0-318 0-648 3:5 0-008 0-853 | 
eS. (0-231 0-705 4-0 0-004 0-856 | 
| 1-4 0-169 0:745 ee 0 0-858 | 


| 
| 


Table 1. The functions J(r) and (7). 





The boundary conditions at an open end can now be derived from 
equation (3) where, in view of the foregoing discussion, all pressures are 
replaced by the corresponding values of the speed of sound, and the 
derivative in the integrand is approximated by a step-wise constant function. 
The integral then breaks down into a series of terms each of which represents 
the contribution at the mth point from all preceding elements of the incident 
wave. One obtains 


n Aa; 
a, = a+ Y —+ [0(r, —7;1)- (7, -7))], (9) 
j=1 77 T§-1 
where the function ®(7) is defined by 
O(r) = | 1(0) d0. (10) 


0 
The function I(r) is known from the previous investigation, and the new 
function ®(7) was derived from it by numerical integration. Values for 
both functions are listed in table 1. ‘These numbers are uncertain within 
a few units in the last decimal place, but this accuracy is entirely adequate 
for wave diagram applications. A plot of ®(7) is also included in figure 1. 
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Equation (9) represents the instantaneous boundary conditions. To 
illustrate the nature of the adjustment process at the duct exit, consider 
an incident compression wave F(r) defined by 

dF(r) ; | 
= (a, —a)/(7, —7) = const. tor 7, <7 <7}, 
- + (11) 


| 





= 0 for tr <7, and 7 > 7}. | 
It is not necessary, in this case, to divide the wave into small elements, since 
the values of Aa;/Az; would be the same for every element. Alternatively, 
one may consider the wave as one element of strength Aa = a,— a, that is 









a t 
‘a ____ INCIDENT WAVE 
_ BOUNDARY 
“ CONDITION 
a : > 





Figure 3. Variation of the flow conditions at the end section of the duct for an 
incident wave consisting of a single element. 

arriving at the exit during the interval Ar = 7,—7). The boundary 

conditions then follow from (9), and are given by 


a= as for tT <7, ] 

| 

a,;— a . 

a = a+ —— O(7-71) for 7) <T< 7, | 
5 lama f (12) 

a,—a f | 

a = a+ —— [O(7-7) -—O(7 -7,)] for 7 > 7}. 

‘ 1 _—s ‘oO J 


These conditions are shown as the heavy line in figure 3. The values 
increase to a maximum given by 
Amax = 4T Seria —T)), (13) 
0 
and then decay asymptotically. If 7, is large, the second term on the 
right side becomes small and amax tends to ad. ‘This is exactly what one 
would expect, for, if the incident wave changes the flow conditions only 
slowly, the boundary conditions do not deviate significantly from their 
steady-flow values. If on the other hand, one considers the limit as 7, — 7 
tends to zero (incident shock wave), one obtains with the aid of equation (10) 
lim ®(7)/7 > (0) = 1, (14) 
7) 


and, therefore, @max = @. 
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The procedure indicated by equation (9), applied to an incident wave 
for which dF/dz is not constant, represents a superposition of the contri- 
butions from all wave elements that have reached the exit since the beginning 
of the wave. ‘The difference of the ®-functions in each term of the sum 
rapidly approaches zero when the smaller argument exceeds about three 
(see figure 2). In practice, one finds that only the most recent three or 
four wave elements contribute significant amounts to the sum. ‘The incre- 
ments after the end of the incident wave are all zero, but the calculations 
must be continued until all the contributions of the earlier elements have 
also decayed to zero. 

This computing procedure is more complicated than the conventional 
one which is based on the simple boundary condition a, = a, but the 
additional time required is not prohibitive. From a, and P,,, one obtains 
u, and O, by means of equations (4), so that the flow conditions and the 
reflected characteristic at the exit are completely determined. 


3. MODIFICATION OF THE PROCEDURE FOR INFLOW 

As discussed in § 1, the procedure derived in §2 applies only to outflow 
from the duct, and some modifications are necessary before it can be applied 
to inflow. ‘These modifications are of an empirical nature, and although 
they seem to be plausible, the results must be checked experimentally 
for a variety of flow conditions. 

First, it is necessary to correct for the variation of the steady-flow 
boundary conditions with the velocity of inflow. ‘The sum on the right side 
of equation (8) represents the deviation of the instantaneous from the 
steady-flow boundary conditions. For inflow, the steady-flow boundary 
conditions are not represented by the constant value a), but by a, ,. ‘The 
latter is related to the steady-flow velocity u, , by the energy equation 
(see Rudinger 1955 a, pp. 63-65) 


Nn, 


2 2 
a ,+u-, = —74@, (1 
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where dp, the speed of sound in the gas surrounding the duct, assumes the 
role of the stagnation condition. Since the Riemann variable P,, of the 
incident characteristic is known, equations (4a) and (15) completely deter- 
mine the flow that would be established if there were no lag in the adjustment 
to the steady-flow boundary conditions (as in the conventional procedure). 
Replacing a, by a,,, in equation (9) ensures, therefore, that the correct 
boundary conditions are obtained when a sufficiently long time has elapsed 
after the end of the incident wave. 

The strength of every wave element Aa; must also be modified to 
account for the’change that the element causes in the steady-flow boundary 
conditions. For an incident expansion wave (to produce inflow) consisting 
of a single element specified by Aa and Ar, the largest deviation from the 
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steady-flow boundary condition occurs at the minimum value of a. "Vhe 
latter is given by the relation, analogous to equation (13), 

Amin = ay 4 — (p(Az), (16) 
where Aa, represents the effective, modified strength of the wave element. 
This relation yields the correct value dmin = a, for large Ar. It must 
also lead to ai, = @+Aa (Aa <0 for an expansion wave) in the 
limit as Av approaches zero. Because of equation (14), one obtains 
then Aa, = Aa—(a,—a,), which shows that the strength of a wave 
element must be decreased by the associated change of the steady-flow 
boundary conditions. For a wave that is formed by a number of 
elements, the effective strength of any one must, therefore, be given by 
Aa;—(a;,—a;,,,). ‘The assumption implied in this analysis is that the 
functions (+) and ®(7), which were derived for the adjustment process 
during outflow, can also be applied for inflow in spite of the different flow 
patterns in the two cases. With these modifications, equation (9) becomes 
Aa; ~ (4j,2 — 41,0) [O(r, —7,4)—(z, —7,)], (17) 


T 
Jj 


a, =Q,.+ > 


| 75 Tj 
which includes the outflow conditions as a special case, since it transforms 
into equation (9) for a; = ap. 

So far only insentropic inflow has been considered, which is permissible 
only if the flow does not separate at the inlet edge of the duct. ‘The losses 
associated with flow separation result in an entropy increase that must be 
taken into account. ‘The characteristics of the incident wave must then 
pass through a region of variable entropy before they reach the end of the 
duct, and the values of P are no longer constant, but vary according to 
(e.g. Rudinger 1955 a, pp. 18, and 39 & 40) 
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where 6. /d7 indicates the derivative in the characteristic direction u+ a, 
and S is a dimensionless entropy (= entropy/(y—1) x specific heat at 
constant pressure) measured from the reference value S, = 0. ‘Therefore 
it is necessary to find how the values of Aa, are affected by the variations 
of P. The effect of entropy gradients can be derived with the aid of figure 4, 
which shows the incident characteristics at the exit points n—1 and n 
and the path of the gas ‘particle’ that separates the regions of zero and 
variable entropies; the ‘particle’ path entering the duct at the point n— 1 
is also indicated. Let P,,_,) and P,,, represent the Riemann variables of 
the two incident characteristics before they enter the region of variable 
entropy. In this region, the strength of the wave element is given by 

Aa, 9 = &.9—- 2n-1,0 = y= nf. * P,-1,): (19) 
The boundary conditions depend only upon waves that actually reach the 
end of the duct and, as stated before, the value of Aa, should therefore 
express the change of the flow conditions that would be produced at the 
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location of the exit if the duct did not terminate there after arrival of the 
preceding wave element. Under these conditions, there would be no 

entropy change between k and » (see figure 4), so that 
P, = P,, and Q), = Q,_.. (20) 
Equations (4) yield 

a, = +y—1)(P,4+9,), anda,_, l(y—1)(P,_,+9,_1), 
and, together with (20), lead to 

Aa, = a@,-4,_, iy — INP. — P,-1): (21) 
Equation (18) can be integrated since a on the right side can be considered 
as a constant. ‘This is always permissible in the cases considered here, 
because changes of @ are small and aS isa sufficiently small quantity, so that 
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Figure +. Reflection of a wave element that first traverses a region of variable entropy 
The dashed lines indicate the path of gas ‘ particles > along which the entropy 


remains constant. 


the errors resulting from the foregoing simplification are negligible. ‘The 
integration yields for the mth characteristic 


i, = | T a, AS, = S,,.0)« (22) 


i] 


If this equation, applied to points k and n—1, is substituted into (21) and 


combined with (5) and (19), one obtains 

Aa, = Aa, 9{1+24(y—1)S,_1}, (23) 
where use has also been made of the condition S; = S,_, (i.e. k and n—1 
lie on the same ‘particle’ path). ‘The intermediate wave-diagram point / 
is determined by P, from equation (22), and QO, from the analogous relation 
O, = O,4+4,-,(S,—S,,_1), where S,; must be found by interpolation. 
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Equation (23) provides the desired correction for the strength of a 
wave element owing to its passage through the region of variable entropy, 
but the correction factor deviates from unity only by an insignificant amount. 
For sonic flow velocity into a sharp-edged duct inlet (Borda mouthpiece), 
the error would be less than 3° for y = 5/3. In all actual cases, 
it would be even smaller, because the inflow velocity must be less than sonic 
for a characteristic to be able to reach the end of the duct. It can therefore 
be concluded that the strength of the incident wave elements may be 
determined from equation (19). 

The values of a,, ,. in equation (17) can be obtained in the conventional 
manner from equations (4a) and (15), but P,, must first be computed 
from P, or P,,y by means of equation (22). ‘This step requires the yet 
unknown value of S,. Some suitable assumptions must be made to deter- 
mine this quantity. The inflowing gas forms a stationary vortex inside 
the duct, followed by a mixing region, after which the flow again fills the 
entire cross section (figure 5). In steady flow into a given inlet configuration, 
the entropy rise is a function of the flow velocity (or its speed of sound, 
since these two variables are related through equation (15)); in nonsteady 
flow, the changes of entropy lag behind those of velocity, because the 
former are a consequence of a mixing process which requires a certain time. 
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Figure 5. Flow into a duct with separation at the inlet edge. 


The following simple assumptions have led to satisfactory agreement 
of the results with experimental observations. ‘The vortex and subsequent 
mixing region are concentrated in the inlet section of the duct, and the 
entropy is given by 


5, @ Baws (24) 

To illustrate the last assumption, consider an incident expansion wave 
consisting of a single element. ‘The speed of sound at the inlet drops 
from its initial value a) to some minimum, 4), after which it gradually rises 
to a,.. ‘The entropy, at the same time, rises from zero to S,, which, 











Reflection of pressure waves 59 


because of the lag, does not reach the level that would correspond to a, 
in steady flow; ultimately, the entropy becomes S,,,._ According to 
equation (24), the assumption is now made that S, = S, .. 

These calculations require the steady-flow relation between entropy 
rise and flow velocity. Sucha relation is readily computed for a sharp-edged 
inlet (Rudinger 1955 a, pp. 71-73), for which it takes the form 

ne is |; ba (40/4 n co) —} i2y y= | (25) 
’ y (49/4, « )? (y-) 
In general, the entropy rise depends on the configuration of the duct, and 
must be given in any particular case by a relation or plot equivalent to 
equation (25). Since a, .. in this relation must be determined from P,,, 
which in turn depends on S,,,, a certain amount of iteration may be 
required for this part of the calculations. 

A further consequence of flow separation at the inlet is the formation 
of a vena contracta (see figure 5). ‘The diameter of the stream tube that 
enters the duct, and in which the flow conditions must adjust themselves 
to the varying boundary conditions, is therefore smaller than the duct 
diameter. Accordingly, the significant length that is used in equation (1) 
to make the time variable nondimensional is assumed to be the diameter of 
the vena contracta D, rather than the duct diameter D. ‘The ratio of the 
cross-sectional area of the vena contracta to the duct area can be computed 
for a sharp-edged inlet as a function of the flow velocity (see, for instance, 
Busemann 1931, p.377), and varies between 0-5 for low velocities (incom- 
pressible flow) and about 0-64 for sonic inflow. Accordingly, D,/D varies 
between about 0:7 and 0-8. Variations of D,/D for other inlet configuations 
might be determined by the method of Jobson (1955), but a constant 
average value D,/D = 0-75 was found to be satisfactory not only for a 
sharp-edged inlet, but also for a flanged inlet with a rounded edge (see § 4). 

The inflow procedures may now be summarized as follows for point 
of the wave diagram. 





(1) Makeanestimatefor S, = S, ,,and compute P,, from equation (22). 


n,% 
The required value of a, 4 or a, is already known from preceding wave 
diagram calculations. 


n,f 

(2) Determine a, ,. from equations (4a) and (15). 

(3) Find S,,.. from a, ,. and the prescribed steady-flow relation between 
entropy and the speed of sound (equation (25) for a sharp-edged inlet or 
its equivalent for other configurations). If S,,,. does not agree with the 
original estimate, repeat steps (1) to (3) with a new value until agreement is 
reached. These three steps also represent the conventional procedures 
which are based on S, = S,,,. and a, = a, ,. 

(4) Determine Aa, from equation (19). 

(5) Compute a, from equation (17), where the value t = agt/0-75D is 
used. (This ‘time variable must be used in equation (17), but it is not 
necessary to use D, as the reference length for the entire wave diagram.) 
The required data for points j < are already known from preceding work. 
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(6) Compute u, and O, from P,, and a, with the aid of equations (4). 
With this step, the flow conditions at point m are completely determined, 
and the reflected characteristic can be plotted by means of the standard 
techniques. 


4, EXPERIMENTAL CHECK OF THE PROCEDURES 
Experimental techniques 

A shock tube constructed of a cold-drawn brass tube of 3-23 in. internal 
diameter was used for the experimental work. ‘The length of the pressurized 
chamber was fixed at 12ft., while the open duct was made up of several 
sections so that its length could be varied between 3 and 15-5 ft. Depending 
on the desired pressures, one or two sheets of photographic film or cello- 
phane were used as diaphragms. 

Two flush-mounted, condenser transducers (Rutishauser model HR3 
with electronic indicator type S'T—12) were available for pressure measure- 
ments. ‘These had a range of 50 1b./in.? and an adequate frequency response 
(the resonance frequency was of the order of 30000c/s). ‘The pressure 
signals were displayed on a dual-beam oscilloscope (Dumont model 322A) 
and photographed on a drum camera (Southern Instruments model M—1020) 
at a film speed of 50 ft./sec. A second oscilloscope provided 0-1 millisecond 
timing marks that were recorded simultaneously with the pressure by means 
of an optical beam splitter arrangement. After each run, a calibration 
signal was recorded on the same film to check the sensitivity of the instru- 
mentation. ‘The calibration of the pressure recording system was linear 
over the range of the experiments. ‘Two switches, built into the camera, 
served to blank the oscilloscope beams, except for one revolution of the 
drum, and to break the shock tube diaphragm at the right instant by means of 
a solenoid-operated needle. 

The wave phenomena in the tube are indicated in figure 2. Comparison 
between theory and experiment is to be made at a point located 1 ft. from 
the open end (I). The incident wave reaches the transducer at A, and the 
reflected wave starts to arrive at B. After this time, the pressure record 
indicates the superposition of thé two waves. <A simple technique to 
determine the pressure distribution of the incident wave alone is to add an 
extension section to the shock tube which delays the arrival of the reflected 
wave sufficiently to record the incident wave over the range of interest in 
a separate experiment. ‘The flow conditions associated with the incident 
wave are completely determined by the pressure record, since a follows 
from the condition of isentropic changes of state (wall friction effects are 
assumed to be small enough to be neglected) 

scans 
i (2) _ (26) 
AQ Po) 
and u follows from the property of the characteristics that O = Qs. 
Equation (4b) yields, therefore, 
u = 2(a—a,)/(y —1). (27) 
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This method can only be used if the experiments are reproducible. 
Otherwise, the incident wave may be found with the aid of a second pressure 
transducer mounted 2 ft. from the open end. At this location (II), the 
incident wave can be recorded directly in the range between A’ and C” 
(see figure 2) which is adequate for the present purpose. 

In the regions A to B, and A’ to B’, both transducers record the incident 
wave alone. Equations (26) and (27) then serve as a check on their relative 
calibrations, since any pressure at location II appears at location | after 
slat t; —t,, = L/(u+a), (28) 
where L is the distance between the pressure transducers. Equation (28) 
applies only to a simple wave for which the characteristics are straight 
lines. ‘his property may also be utilized to detect condensation of water 
vapour in the region between the transducers, since one would no longer 
have a simple wave (see § 5). 

Once the incident wave is determined by means of one of the two methods 
described, its reflection from the open end may be computed on the basis 
of the conventional and the new procedures. ‘lhe interaction of these 
waves with the incident wave is then determined at location I, and the 
results are compared with the pressure record at this location. 


Reflection of a compression wave 

‘he compression wave was produced by placing a loosely fitting light 
piston into the shock tube a few inches downstream of the diaphragm ; 
when the latter broke, the air pressure accelerated the piston thereby 
producing the compression wave. Proper selection of the weight of the 
piston, the pressure of the driving air, and the length of the duct, ensured 
that the compression wave had sufhcient strength without steepening to 
a shock wave inside the duct, and that the wave reflected from the open end 
passed the points of pressure measurement before the piston arrived there. 
The piston was constructed of two { in. discs of masonite separated by 
five wooden spacing rods of 32 in. length; its weight was about { lb. 
The initial pressure difference across the shock tube diaphragm was 
50 lb./in.?, and the open end of the tube was 15-5 ft. from the diaphragm. 
A new piston had to be used for every run, and instead of relying on the 
reproducibility of the experiments, the incident wave was determined by 
means of the described technique based on two pressure transducers. 
The results obtained in this manner are presented in figure 6. Although 
the incident wave has just steepened sufficiently to form a weak leading 
shock wave, most of the pressure rise in the incident wave takes place 
gradually. ‘The presence of a shock front requires that the term 
Asiocx (7 —7)) (from (8)) be added to the boundary conditions given 
by (9). The experimental results for the reflected wave are in good 
agreement with the data based on the new procedure, while the conventional 
technique indicates too rapid a pressure drop. 
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Reflection of an expansion wave 

The shock-tube chamber was evacuated instead of pressurized in these 
experiments, and to obtain sufficiently steep expansion waves, the distance 
from the diaphragm to the open end was made rather short; it was 3-25 
or 4:25 ft. depending on whether one or two pressure transducers were 
used. Initial pressure differences across the diaphragm were varied between 
—2 and — 13 |b./in.? 
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Figure 6. Reflection of a compression wave. Comparison of theory and experiment 


at 1 ft. from the open end of a duct of 3:23 in. diameter. ‘The time is measured 
from the arrival of the incident wave at 2 ft. from the end. Experimental data : 
A incident wave; 0 interaction of incident and reflected waves. Computed 
data: @ new boundary conditions, [j steady-flow boundary conditions. The 
computed ‘ effective ’ exit pressure ratio is also plotted. 


A sufficiently weak incident wave passes the pressure transducer 
completely before the reflected wave arrives so that both waves can be 
measured directly in one experiment with a single transducer. For stronger 
waves, the two waves interact at the location of the transducer and the 
incident wave must be found by other means. Since the pressure records 
were reproducible, it was possible to use the method of the tube extension 
described in the foregoing, although the alternative method based on two 
pressure transducers was also employed. 

All experiments, with one exception, were carried out with a sharp-edged 
inlet (Borda mouthpiece) because the entropy rise is then readily determined. 
The results for this configuration and for initial diaphragm pressure 


differences of —2, —5, and — 13 lb./in.? are shown in figure 7 as curves a, 
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b and c¢ respectively. It is seen that the customarily used steady-flow 
boundary conditions predict a pressure rise in the reflected wave that is 
faster than the one actually observed, while the new procedure leads to 
good agreement with the experimental data for the two weaker waves. 
For the strongest wave c, the agreement is still good at the beginning of 
the reflected wave, but the computed values gradually rise faster than the 
experimental cnes (see § 5). 
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Figure 7. Reflection of expansion waves. Comparison of theory and experiment 
at 1 ft. from the open end of a duct of 3-23 in. diameter; a, 6, c sharp-edged 
inlet for initial diaphragm pressures of —2, —5, and —13 lb./in.? respectively; 
d flanged inlet with rounded edge; initial diaphragm pressure —5 lb./in.2 (See 
figure 6 for explanation of the symbols.) 


One experiment was carried out with an inlet configuration consisting 
of a flanged duct with a rounded inlet edge, and an initial diaphragm pressure 
difference of —5 Ib./in.?. From the measured steady-flow pressures after 
the incident and reflected waves, respectively, the entropy rise for the final 
steady-flow velocity could be computed and was found to be 61% of that of a 
sharp-edged inlet for the same flow velocity. ‘The entropy rise at the other 
flow velocities was then assumed to be given by the same fraction of the 
corresponding value for a sharp-edged inlet. The results for this experiment, 
plotted as curve d in figure 7, show the same good agreement between the 
computed values based on the new boundary conditions and the experimental 
data. (The incident waves for cases 6 and dare different in spite of the same 
initial pressure conditions, because the diaphragm materials were not the 
same in these experiments ; photographic film was used for d and cellophane 
for d. ‘These variations do not affect the results, since the evaluation of the 
experiments is always based on the individual incident wave.) 


: 5. DiIscussIoN AND CONCLUSIONS 


The previously derived boundary conditions for the reflection of shock 
waves from an open end have been extended to other waves. As long as 
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subsonic outflow from the duct is maintained no further assumptions are 
required. Since the analysis for shock waves led to good agreement with 
experimental observations, one should expect similar agreement for other 
waves for which the disturbance at the duct exit is less violent. ‘This is 
indeed observed, as indicated by figure 6. As additional interesting 
information, the variations of the effective exit pressure were obtained 
from the wave diagram that had to be prepared for the evaluation of the 
experiment, and these are also plotted in figure 6. Itis seen that this pressure, 
after an appreciable disturbance, tends to readjust itself to the steady-flow 
value pp», but levels off slightly above py where it remains almost constant 
as long as the incident wave ‘tries’ to raise the pressure. 

The computing procedures had to be modified to make them applicable 
to inflow. Because of the empirical nature of the modifications, it was 
important to verify their consequences for a variety of flow conditions. 
The described modifications are the ones found to be most satisfactory, 
and the good agreements between calculated and experimental data to 
which they lead is demonstrated in figure 7. A disagreement was found 
only for the strongest waves (c), and then only for the later portion of the 
reflected wave when the observed pressures rise higher than those computed 
on the basis of the steady-flow boundary conditions. Since the difference 
between the two theoretical curves is caused by the lag in the establishment 
of the steady-flow boundary conditions, one must conclude that the dis- 
crepancy cannot be caused by an error in the evaluation of this lag. Any 
uncertainty about the incident wave that might have been caused by the 
effects of condensation of water vapour in the duct was also ruled out by the 
described technique of using two pressure transducers. ‘lhe recorded 
pressures were found to be related to each other exactly in the manner 
required by equations (26) to (28) for a simple wave, which indicated the 
absence of condensation effects. ‘This finding seems to leave wall friction 
effects as the only cause for the observed discrepancy between theory and 
experiment. In this case, the flow in the vena contracta has already become 
sonic, and the Mach number of the flow exceeds 0-5 further inside the duct. 
A crude analysis of the friction effects by a linearization method, similar 
to that of T'rimpi & Cohen (1955), indicates that a friction coefficient several 
times that of its steady-flow value could indeed explain the discrepancy. 
This result appears reasonable in view of the uncertainty of the friction 
coefficient for nonsteady flow in the vicinity of a sharp-edged inlet. 

It is apparent from figures 6 and 7 that the actual flow conditions lag 
behind those computed in the conventional manner by about the time in 
which a sound wave travels one or two duct diameters. ‘The corresponding 
lag time for shock reflection is (see figure 1) the travel time for about three 
duct diameters. 

An interesting consequence of the lag is the elimination of certain 
discontinuities of the incident wave from the reflected wave. Figures 6 
and 7, for instance, show clearly that the derivative dp/dr has a discontinuity 
at the head of the incident wave which reappears at the head of the reflected 
wave only when the conventional boundary conditions are used, while the 
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George Rudinger, The reflection of pressure waves of 
finite amplitude from an open end of a duct, Plate |. 
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Figure 8. Schlieren streak record of the reflection of an expansion wave from an open end of a 
3x3 in. shock tube. Note the apparent absence of a reflected wave. (Courtesy of Dr I. I. Glass, 
University of Toronto.) 
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new boundary conditions, in agreement with the experimental data, indicate 
a gradual change of the slope of the curves. A general proof for this 
observation can be obtained by differentiating equation (3) with respect 
to 7 which yields 
dp rt)  dF\(r) ‘7 dkF(@) dl(z—8@) 
lz dr ss lle —d0-0 dr 
ri ( spies. a7 
where /(0) = 1. ‘The integral in this relation is a continuous function of 7 
if F(z) is continuous, so that any discontinuities of dp,/dz are fully accounted 
for by those of the derivative of the incident wave alone. Consequently, 
no discontinuities appear in the reflected wave. ‘The significance of this 
conclusion becomes evident if one considers that a schlieren photograph 


dd, (29) 





of a gas flow does not record gas densities but density gradients. ‘lhe 
discontinuity of the gradient at the head of an expansion wave in a shock 
tube is therefore readily photographed, while the head of the wave reflected 
from an open end, having no such discontinuity, would appear on the record 
only faintly, if at all. ‘This phenomenon was actually observed at the 
University of ‘Toronto (Glass & Patterson 1955), where schlieren streak 
records of this reflection process were obtained. ‘The relevant portion 
of such a record is reproduced in figure 8 (plate 1), and the originally 
puzzling observation of an apparently missing reflected wave is now quite 
understandable. 

In conclusion, it may be stated that new boundary conditions for wave 
reflection from an open end have been derived which lead to a better agree- 
ment with experimental observations than the customarily used steady-flow 
boundary conditions. ‘The new procedures take into account that the 
steady-flow boundary conditions are not instantaneously established, and 
that the adjustment process is continually modified by the incident waves. 
The actual flow conditions lag behind those computed in the conventional 
manner and, although the lag times are small, they may occasionally become 
significant. ‘Che new computing procedures are more complicated than 
the conventional ones, but the additional time required for their application 
is not prohibitive. On the other hand, one must also consider that the 
resulting refinement may not be warranted because of other errors introduced 
into a wave diagram by some simplifying assumptions, such as the neglect 
of wall friction at high flow velocities. It has also been shown that certain 
discontinuities of the incident wave are eliminated from the reflected wave. 
This improved understanding of the reflection phenomena made it possible 
to explain certain, previously puzzling, experimental observations. 


‘l'his work was sponsored by Project Squid which is supported by the 
Office of Naval Research under Contract N6o-ori-105, T.O.1II]. ‘The 
author is greatly indebted to Mrs Angela Chang who prepared the many 
wave diagrams that were needed to test various modifications of the 
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fications, and to Mr L. M. Somers who carried out the shock-tube 
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SUMMARY 

This paper is concerned with some statistical properties of the 
displacement of a marked fluid particle released from a given 
position in a turbulent shear flow, and in particular with the 
dispersion about the mean position after along time. It is known 
that the dispersion takes a simple asymptotic form when the particle 
velocity is a stationary random function of time, and that analogous 
results are obtainable when the particle velocity can be transformed 
to a stationary random function by suitable stretching of the 
velocity and time scales. ‘The basic hypothesis of the paper is that, 
in steady free turbulent shear flows which are generated at a point 
and have a similar structure at different stations downstream, the 
velocity of a fluid particle exhibits a corresponding Lagrangian 
similarity and can be so transformed toa stationary random function. 

The velocity and time scales characterizing the motion of a 
fluid particle at time ¢ after release at the origin are determined 
in terms of the powers with which the Eulerian length and velocity 
scales of the turbulence vary with distance w from the origin. ‘The 
time scale has the same dependence on ¢ for all jets, wakes and 
mixing layers (and also for decaying homogeneous turbulence) 
possessing the usual kind of Eulerian similarity. ‘The dispersion 
of a particle in the longitudinal or mean-flow direction (and likewise 
that in the lateral direction in cases of two-dimensional mean flow) 
is found to vary with ¢ in such a way as to be proportional to the 
thickness of the shear layer at the mean position of the particle. 
The way in which the maximum value of the mean concentration 
of marked fluid falls off with ¢ (for release of a single particle) or 


with x (for continuous release) is also found. 


1. INTRODUCTION 

The purpose of this paper is to show that the well-known analysis of 
the dispersion of a marked fluid particle due to turbulence can be adapted 
to a discussion of diffusion in (statistically) steady turbulent jets, wakes and 
mixing layers. ‘These cases, which hitherto have lain outside the scope of 
the available analysis, are distinguished by their property that the mean 
velocity is only approximately unidirectional. As a preliminary to an 
explanation of the modification required for these new cases, the existing 
analysis will be described briefly. 
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‘The object, in either a theoretical or an experimental investigation of 
turbulent diffusion of some conserved quantity not subject to molecular 
diffusion, is essentially to determine the statistical history of a volume of 
marked fluid which occupies a given position at an initial instant. If only 
the probability of finding marked fluid at any point and at any subsequent 
time (which is effectively the local mean concentration of the quantity 
undergoing diffusion) is required, it is sufficient to consider each volume 
element of the marked fluid in isolation. Only this limited aspect of 
turbulent diffusion will be considered herein. We are thus concerned with 
the statistical history of a single small element of fluid whose position is 
given at an initial instant. A description of the place of such ‘ one-particle 
analysis’ in the general framework of the theory of turbulent diffusion will 
be found in a recent survey (Batchelor & ‘Townsend 1956). 

It will be enough for us to consider one component of the vector position 
of the fluid particle or volume element under discussion. ‘Thus, we wish 
to determine the probability distribution of the displacement X(t) which the 
particle has undergone during the time interval ¢ subsequent to the initial 
instant at which its position was given. ‘The displacement X(t) is a random 
quantity, taking different values for each trial or realization of the turbulent 
How, and averages will be regarded as referring to an ensemble of such 
realizations. From a practical point of view, the most useful parameters 
of the probability distribution of X(t) are the mean displacement 

t 


y= | u(t’) dt’, 
0 
where u(t) is (one component of) the particle velocity at time ¢, and the 
dispersion D,(t) about the mean displacement, where 
[D,(t)}? = [X(t)— X(OP. 
G. |. ‘Taylor (1921) pointed out that there were advantages in thinking of 
the dispersion in terms of the particle velocity by means of the relation 


= = 2[u(t) —u(t)][X()- X(0)] 


=2 | [u(t)—u(t)][u(t’)—u(t’)] dt’. 


/ 0 











‘L'aylor also noted—and this is the vital point in this well-known analysis— 
that it was possible to obtain some simple and definite results in the case in 
which u(t) is a stationary random function of ¢, i.e. when the statistical 
properties of the particle velocity are constant in time. We then have 





X(t) = tu, 
[u(t) — u]{u(t’)—u] = R(t—-t’), 
dD- - : 
and = 2 a R(&) dé, 


R being the covariance function of the fluctuation in particle velocity. 
It will normally happen, in cases of well defined turbulent flow, that R(é) + 0 
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as € -> 00 and that at least the first few integral moments'of R(é) converge. 
Consequently we have the important result that, as t > «, 
D? + 2At—2B, 


where (as found from integration by parts), 


A= R@dé, B= { ERG a. 
) 


“OU J { 
Moreover, since the displacement can be written in the form 


tle pra 


X(t)= > | u(t’)dt’, 


n=1 / (n—1l)a 

which is a series of random terms of equal variance, it follows from a plausible 
application (although a non-rigorous one, because the terms of the series 
are not statistically independent) of the central limit theorem that, as 
t-—> © (with « remaining constant), the probability distribution of X(t) 
tends to the normal or Gaussian form (Batchelor 1949). When X(t) is 
normally distributed, the probability of finding the fluid particle within 
a given small range of positions at time ¢ satisfies a Fickian or heat conduction 
type of equationand the above integral denoted by A can then be interpreted as 
a diffusivity, or coefficient of diffusion, with dimensions (velocity) x (length). 

In short, it is possible to obtain some useful results, and in particular 
to determine the asymptotic dependence of the dispersion on the time, 
when the random velocity of the particle is a stationary function of time. 
We are thus led to look about for turbulent flows in which the particle 
velocity has this simple property—or, failing stationarity, some related 
property which might allow a modified form of the above analysis to be 


applied. 


2. [WO PREVIOUS APPLICATIONS OF ONE-PARTICLE ANALYSIS 

The first case of turbulent motion in which a study of particle dispersion 
was made was decaying turbulence generated by a grid placed in a uniform 
stream (Taylor 1935), and this remained the only one for many years. 
As a consequence of the decay of the turbulence, the velocity of a fluid 
particle is not here a stationary random function of time (not even 
approximately), and some modification of the above analysis is necessary 
before it can be applied. Provided it may be assumed that the turbulence 
preserves its structure as it decays (so that all functions characterizing the 
turbulence have the same form, although possibly with different length, 
time and velocity scales, at all relevant stages of the decay), a suitable 
modification of the diffusion analysis can be made to allow for the decay 
(Batchelor 1952; Batchelor & Townsend 1956). The essential point of 
the modification is to make a mathematical transformation of the particle 
velocity —expanding the velocity scale and contracting the time scale more 
and more as the decay proceeds—such that it becomes a stationary random 
function of a new variable related to time ina known way. (Taylor (1935) 
attempted to allow for the decay of the turbulence, but his method is not 
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wholly correct inasmuch as it takes into account only the change in the 
velocity scale and ignores the possibility of an independent change in the 
time scale of the motion.) ‘The modification that will be adopted to make 
the diffusion analysis applicable to free turbulent shear flows is similar in 
principle to that used for decaying turbulence behind a grid. 

A few years ago, when Sir Geoffrey Taylor was conducting his 
experiments on the longitudinal extension of a finite volume of salt solution 
injected into the turbulent flow of water along a circular pipe (‘Taylor 1954), 
I realized that turbulent flow along a cylindrical pipe is a case to which the 
above one-particle analysis can be applied without the need for any modifica- 
tion. ‘The two properties of the turbulent flow in a pipe that together 
ensure that the velocity of a particle is a stationary random function of time 
are: (a) the fluid particle is constrained by the walls always to lie within the 
turbulent flow inside the pipe; (4) the turbulence has the same statistical 
properties at all cross-sections of the pipe. Even though the fluid particle 
may move into different parts of the pipe cross-section, including the slow 
moving layer near the wall, the statistical properties of the particle velocity 
do not change as it moves along the pipe, and the above analysis of diffusion 
(with w and wu referring to the longitudinal or axial direction) is immediately 
applicable. Some observations of the times of travel of small solid spheres 
between two distant cross-sections of a circular pipe have been analysed with 
the aid of this diffusion analysis (Batchelor, Binnie & Phillips 1955). 

We are now in a position to consider the modification necessary to render 
the diffusion analysis applicable to the class of free turbulent shear flows, 
of which jets, wakes, and mixing layers are examples. Free turbulent shear 
Hows differ from turbulent flow in a pipe, firstly in that the turbulence is 
not bounded by a rigid wall but adjoins non-turbulent fluid of uniform mean 
velocity, and secondly, in that, as a consequence of the absence of the 
rigid cylindrical wall, the turbulence changes from one cross-section to 
another, the simplest manifestation of this change being an increase in the 
width of the turbulent flow with increase of distance along the streamlines 
of the mean flow. So far as this second difference is concerned, it will 
often happen, in cases of free turbulent shear flow, that the turbulence 
preserves its structure as it changes with distance downstream, and that the 
change in the turbulence is confined to changes in the length, time and 
velocity scales of the motion; such a state of self-preservation is usually 
set up at a sufficient distance from the source or origin of the turbulence*. 

* Turbulent boundary layers on rigid walls have one ‘ free’ boundary and thereby 
might be thought to qualify for inclusion within the group of shear flows under dis- 
cussion. However, boundary layers on uniformly rough walls are different from the 
other examples named inasmuch as a self-preserving state of the turbulence is attained 
here only in a certain limited approximate sense, even when the Reynolds number is 
large. (In the rather artificial case of a uniform stream passing over a rigid wall with 
linearly increasing roughness height, exact self-preservation of the boundary layer 
turbulence is possible; the diffusion laws are then the same as those to be described 
later for the case of a mixing layer with zero velocity in the stream on one side.) In 
what follows, only those free turbulent shear flows which attain a self-preserving state 
exactly (possible only asymptotically, i.e. far downstream) will be considered. 
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‘ We have here a situation which is a little like that in the case of diffusion 

‘J in decaying turbulence behind a grid, and a natural suggestion is to employ 

“ the same trick of transforming the velocity of a fluid particle in such a way 

™ that it becomes a stationary random function. However, before exploring 
this suggestion, there is a question arising out of the first of the above 

» differences between free turbulent shear flow and confined turbulent shear 

. flow which must be answered. 

) 

e 3. CAN A FLUID PARTICLE ESCAPE FROM THE REGION OF TURBULENT 


; MOTION ? 
j As remarked above, the presence of the walls of a pipe ensures that a 
fluid particle remains within the region of turbulent flow. Is there a 
similar guarantee, in the case of free turbulent shear flows, that a particle 
remains always in the region of turbulent flow? Or is it possible that some 
fluid elements are ejected from the region of turbulent flow into the 
surrounding non-turbulent fluid and subsequently take no part in the 
turbulent motion? It is necessary to settle this question before proceeding 
to an analysis of turbulent diffusion, because, if the possibility of such an 
escape from a region of turbulent flow to a region of non-turbulent flow should 
exist, it would have a profound effect on the stat:stical properties of the 
particle velocity. ‘lhe velocity of the fluid particle would fall rapidly to 
zero, and might stay at that value since a return to the region of turbulent 
flow would not be inevitable. (‘This is in contrast to the case of flow in 
a pipe; the viscous sub-layer near the walls may perhaps be regarded here 
as a region of non-turbulent flow in the particular sense that inertia forces 
are not important, but the velocity of a particle remains random in the 
viscous sub-layer and, in view of the geometrical contraint of the pipe wall, 
there is statistical certainty that a particle in the viscous sub-layer will 
eventually return to the central region of the pipe.) In the event of the 
probability of escape from the region of turbulent flow being finite, it would 
be impossible to transform the particle velocity into a stationary random 
function by simple adjustment of the velocity and time scales. 

Fortunately, the available evidence about the nature of free turbulent 
shear flows suggests that such escape does not occur (see Corrsin & Kistler 
1954). It has been known for some years that free turbulent shear flows 
are characterized, instantaneously, by a sharp boundary, of irregular and 
random form, separating a central region of turbulent motion from an outer 
region in which the flow is non-turbulent in the sense that the vorticity is 
zero. Something of the mechanism by which this instantaneous boundary 
remains sharp against the smearing action of viscosity is known. It seems 
that vorticity is diffused, by the action of viscosity, from the central to the 
outer region, and that the high rate of stretching of vortex lines in the central 
region rapidly increases the magnitude of the vorticity to some high equili- 
brium level as,soon it is made finite by viscous diffusion; in this way the 
sharp boundary propagates relative to the non-turbulent fluid in the outer 
region. What is important in the present connection is that the boundary 
always advances into the non-turbulent fluid, thus acting as a valve which 
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allows non-turbulent fluid to pass into the central region by mixing or 
entrainment and to be converted to turbulent fluid, but which does not allow 
fluid to pass out of the central region. 

It seems, therefore, that, if a fluid element is once inside the central 
region of turbulent motion, it remains within it. ‘The general way in which 
the velocity of a fluid particle in the central region changes with time will 
thus be related to the fact that, as the particle moves downstream (always 
remaining within the central region), it is subject to the influence of turbulent 
motion whose length and time scales are continually changing. We now 
consider how the velocity of a fluid particle in the central region of turbulent 
motion may be transformed into a stationary random function. 


4, 'TRANSFORMATION OF THE PARTICLE VELOCITY FOR FLOWS WITH 
SIMILARITY 

We shall suppose that each of the turbulent shear flows concerned is 
steady and has the same structure at different distances x downstream from 
some virtual origin. ‘The only change in the statistical properties of the 
flow (including the variation of the Eulerian mean velocity, although not 
necessarily its absolute magnitude) at different stations downstream is 
a change of the length (L) and velocity (V) scales of the motion. It will 
also be assumed that these length and velocity scales are proportional to 
powers of x, Le. 

L(x) oc x?, iz) 24, (4.1) 
the corresponding variation of the time scale of the motion being as x?*%. 
Both these assumptions are usually valid at sufficiently large Reynolds 
numbers and at sufficiently large distances from the real origin for the 
common cases of free turbulent shear flow which are not influenced by 
rigid boundaries. 

An element of fluid is carried downstream, and the statistical properties 
of its velocity u(t) change with ¢ as a consequence of the variation of 
the properties of the Eulerian velocity field with position x. Whatever 
may be the distance of the element downstream, it finds itself surrounded 
always by turbulence of the same statistical form. Thus the fluctuations 
of the velocity of the fluid particle have a (statistically) similar form 
at different times. We are therefore led to make the hypothesis that 
the particle velocity can be transformed to a stationary random function 
by suitable adjustment of the velocity and time scales. (More precisely, 
it is not u(t), but u(t)— Uy, which is so transformable, where U, is a velocity 
of translation of the whole field which has no effect on the turbulence and 
serves only to give a frame of reference with respect to which the flow is 
statistically steady; Ug is finite in the case of a wake, and is usually zero 
otherwise.) Analytically, the hypothesis is that 


[u(t) — Up]/w(t) = F(n), (4.2) 


where F(7) is a stationary random function of a new variable 7 given by 
dyn & dt/7(t), (4.3) 
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and w(t) and z(t) are the velocity and time scales of the statistical properties 
of the particle’s motion at time ¢. ‘The origin of t can be taken as the instant 
of release of the particle from a fixed position, and it will be necessary for 
validity of the hypothesis that t be large enough for the exact circumstances 
of the release to be ‘forgotten’. ‘The hypothesis represented by (4.2) 
and (4.3) is no more than a similarity hypothesis for the Lagrangian features 
of the turbulence, exactly analogous to the better known similarity hypothesis 
for the Eulerian features. It is probable that the one similarity hypothesis 
is a strict consequence of the other, but attempts to prove this encounter 
the usual difficulties of relating Lagrangian and Eulerian features of flow. 

We have now to determine w(t) and 7(t), making use of the information 
that the velocity and time scales of the Eulerian properties of the turbulence 
at distance x downstream are proportional to x“ and x’*¢ respectively. 
Now the mean distance downstream of the particle (which will be supposed, 
for the sake of simplicity of algebraic form of the relations that follow, to 
have been released at x = 0, i.e. at the origin of the turbulence) at time ¢ 
is X(t), and the velocity and time scales of the Eulerian features of the 
turbulence at the position x = X(t) will presumably be also the velocity 
and time scales of the statistical properties of the particle motion at time f. 
Thus we have a _ 

w(t) oc [X(t)}-4%, (t) oc [X(t) PP *4. (4.4) 
The mean velocity of the particle is itself one of the statistical quantities 
included in the similarity hypothesis, and the dependence of X(t) on t can 
be determined from __ 
XO) _ Gy = Uy+ Fuld, (4.5) 
dt 
in which F(») is a constant. 

Equations (4.2) and (4.3), together with the auxiliary relations (4.4) 
and (4.5), define the transformation that enables a modified form of the 
diffusion analysis described earlier to be applied to free turbulent flows 
whose properties change with distance downstream. 

The form of X(t), as determined by integration of (4.5), depends on 
whether U, is zero or not. It is convenient therefore to give separate 
consideration to these two cases in the next two sections. 


5. DIFFUSION IN JET-TYPE FLOWS 
We include under this heading all those steady free turbulent shear flows 
in which, as in the typical case of a jet discharging into stationary fluid, the 
absolute value of the velocity at any point conforms to the similarity laws 
represented by (4.1), so that Uy = 0. Integration of (4.5), with the aid of 
(4.4), then gives 
X(t) oc PG+9, (5.1) 
the constant of integration being determined by the condition that the fluid 
particle is released at the position x = 0. 
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The new independent variable 7 is thus defined by 
dyn = (P+ +9 dt, 
the arbitrary multiplicative constant being chosen as unity for convenience. 
Integration leads to two different forms of relation between 7 and ¢ according 
asp =lorp-+1. It is not necessary to consider both cases, because, for 
all the types of flow concerned in this section, similarity of structure of the 
flow at different distances downstream is possible only with p= 1. This 
may be seen from a comparison of the two terms U 0U/ox and euv/oy which 
occur in the Eulerian equations of mean motion (the symbols having their 
usual meanings and not those applicable elsewhere in this paper). When 
both U? and uv have the form required for (Eulerian) similarity, viz. 
v *4 x function of y/x”, 
the equation can be satisfied only if p= 1. Jets and mixing layers (with 
zero velocity on one side) are thus straight-sided, the turbulence being 
contained in either a conical or a wedge-shaped region. 
The relation between 7» and ¢ is thus 
y = logt. (5.2) 
The fluctuation of the particle velocity about its mean value is 
u(t) —u(t) oc t-4+9[ F(n)F —(n)] 
=f “as YF(n), Say, 
where f(7) is likewise a stationary random function of y, and the fluctuation 
of the particle position about its mean value is 
” t 
X(t)— X(t) = | U-UG+DF(y’) dt’ 
“0 


log t 
= | evit+nf(y’) dy’. (5.3) 
We note in passing that the presence of the weighting factor in this integral 
will lead to the value of X(t)— X(t) being dominated by a finite portion of 
the range of integration, even when t-> 0; in these circumstances, it 
would not follow from the central limit theorem that the probability distribu- 
tion of X(t)— X(t) tends to the Gaussian form as t > ©. 
The relation from which the statistical dispersion of the particle is 
determined now becomes 
dD? ‘logt 


“o 


— f0-M/A+a | en +MR(L) dC, ( 


wa 
> 
~- 


where ¢ = »—7’ and R(C) is the covariance function of f(y). If the position 
of release of the particle had been taken as some finite value of x, the upper 
terminal of this integral would have been finite; however, the integral is 
convergent and the asymptotic form of the dispersion (as t > ©) would 








“I 
ui 
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still be given by the above relation. Integration of (5.4) gives 

D,(t) « t+ oc X(t). (5.5) 
The striking feature of this result is that the longitudinal dispersion of the 
marked fluid particle increases in proportion to the thickness of the shear 
layer at the mean position of the particle. Comment on this feature is 
postponed until § 8. 

If it happens that the probability distribution of the displacement of the 
fluid particle has a Gaussian form (which may be so in some cases, although 
there are as yet insufficient grounds for expecting the Gaussian form in 
general), the quantity $}dD=/dt can be interpreted as a diffusivity. We 
might then account for the power of ¢ in the expression for }dD*/dt 
by noting that ¢!"*+® is a measure of the mean distance downstream 
of the fluid particle at time ¢, and that representative length and velocity 
scales of the turbulent diffusing motions at distance x downstream are 
proportional to x and x~¢ respectively. 


6. DIFFUSION IN WAKE-TYPE FLOWS 
In the case of a wake behind a body held fixed in an otherwise uniform 
stream of speed U,, only the velocities in the shear flow relative to U, 
conform to the self-preservation laws, and the first term on the right hand 
side of (4.5) is non-zero. Moreover, a turbulent wake has a self- preserving 
structure only when the variations of velocity across the wake are small 
compared with U,, so that we might as well neglect the second term on the 
right hand side of this equation. With this approximation that the mean 
position of the particle is the same as if it travelled with the speed of the 
undisturbed stream, i.e. that _ 
Ald) = Ugh, (6.1) 
we have dy =t?-“dt 
(the multiplicative constant again being put equal to unity). Again integra- 
tion gives two different forms of the relation between » and ¢ according as 
p+q=1orp+q #1, and again an additional condition for self-preservation 
of the turbulence allows only one of these possibilities. As before, we 
compare the ways in which the terms U 0U/dx and éuv/dy in the Eulerian 
equation of motion depend on x; on noting that the appropriate approximate 
form of UdU/dx is U, 0U/ex in the present case, we find 
pt+q=1. 
Thus the relation between » and ¢ is again 
n = logt. 
The fluctuation of the particle velocity about its mean value is now 


u(t) —u(t) oc (Uyt) “[F() — F(y)) 
=t%f(m), say, 
and the fluctuation of the particle position about its mean value is 


-logt 


X(t)— X(t) = | f(y!) dry. (6.2) 
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The relation from which the dispersion can be determined is 





dD: togt 
a eee | eT hr) an 
= pa [eto R(DAL, (6.3) 
showing that D,(t) oc 8-4 o LX(é)). (6.4) 


The general remarks in the preceding section are also relevant here. 


7. LATERAL DIFFUSION IN TWO-DIMENSIONAL FREE TURBULENT SHEAR FLOWS 
The notation and wording of the preceding sections have been chosen 
to refer to diffusion in the longitudinal direction, i.e. in the direction of the 
mean flow. However it is evident that, in cases in which the shear flow is 
statistically two-dimensional, the above analysis can also be regarded as 
describing diffusion in the lateral y-direction (that in which the turbulence 
is stationary). ‘The convection of the marked particle in the longitudinal 
direction by the mean flow again determines the way in which the properties 
of the turbulence in the neighbourhood of the particle change with time, 
and lateral diffusion follows laws of the same formas for longitudinal diffusion, 
the only difference being that there is no mean displacement of the particle 
in the lateral direction. ‘Thus, the lateral component of the velocity of the 
particle is of the form 
(t-%+%o(n)_ for jet-type flows, 
a t~%9(7) for wake-type flows, 


where g is a stationary random function of y (which is related to ¢ as in the 
two preceding sections) with zeromean. Then, if Y(t) is the lateral displace- 
ment of the particle from its initial position in time ¢, the mean value of Y(t) 
is zero, and the dispersion about the mean is given asymptotically, as t > ©, 
by 

{ -o 
| 4d—q/+q) | e7e/ +9 S(C) dé for jet-type flows, 
dD: | 0 

XL 


dt 7 “o. 
th-2a | eh S(L) dl for wake-type flows, 


where S(C) is the covariance function of g(y). ‘Thus, in all cases, 
D,(t) « D,(t). (7.1) 


‘The proportionality of the longitudinal and lateral dispersions is a simple 
consequence of the fact that the transformations needed to convert the 
longitudinal and lateral components of the particle velocity into stationary 
random functions are of the same form. However it should be kept in 
mind that, as a result of the shearing action associated with the distribution 
of mean Eulerian velocity, the fluctuations in u(t) are greater in magnitude, 
and possibly persist for a longer time, than those in v(t), so that D_,(t) may 


Diffusion in free turbulent shear flows 77 


be a good deal larger numerically than D,(t). ‘The increase in thickness of 
the shear layer with distance downstream is essentially a process of diffusive 
spreading of fluid marked with finite vorticity and, since the two lateral 
components of velocity have comparable intensities in two-dimensional 
shear layers, D,(t) is likely to be of the same general magnitude as the 
thickness of the shear layer at x = X(t). 

The lateral spread of positions of the marked particle would be of 
particular interest in a case in which marked particles are released at the 
origin at regular intervals, one after the other, or when the release is con- 
tinuous, as it might be when dissolved salt is used as the method of marking 
the fluid. ‘The boundary, in the (x, y)-plane, of the ‘ wake’ of marked fluid, 
defined as the curve on which the mean density of marked fluid is some low 
arbitrarily chosen fraction of the mean density at y = 0 for the same value 
of x, is then a curve whose ordinate is proportional to D,(t) when the abscissa 
has the value X(t). ‘Thus the marked fluid is bounded by the curve 


y ox=x? for jet-type flows, 
or y cxlt*= x? for wake-type flows. 


The thickness of the ‘ wake’ of marked fluid in the y-direction is proportional, 
in all cases, to the thickness of the turbulent shear layer in the z-direction. 
The marked fluid spreads out in the z-direction as rapidly as the growth in 
thickness of the shear layer allows it, so that when marked fluid is released 
continuously the ‘ wake’ of marked fluid has the same cross-sectional shape 
at all values of x. 


8. SOME COMMON FEATURES OF THE ABOVE RESULTS 

Despite the apparent differences between the various cases of free 

turbulent shear flow, some aspects of the results obtained in §§5-7 are 

common to them all. ‘Two of these common features in particular call 

for notice in view of their fundamental character. ‘The first is the relation 

between ¢ and the new variable 7. For both jet-type and wake-type flows, 
this relation was found to be 

y= log t; (8.1) 


this was also the relation found in the case of decaying turbulence behind 
a grid (Batchelor & Townsend 1956). Going back one step, the common 
feature of all these cases of turbulent flow which leads to this logarithmic 
relation is the fact that the time scale characterizing the motion of a fluid 
particle is proportional to t, where ¢ is measured from the instant at which 
the particle is released (at the point at which the turbulence originates). 

It is remarkable that in such widely different types of turbulent flow 
as jets, mixed layers, wakes, and decaying grid turbulence, the time scale 
of oscillations in the velocity of a fluid particle should always increase 
linearly with ¢. ‘The distinctive property common to all these developing 
or decaying turbulent flows is their self-preservation, or similarity of 
structure at different stations normal to the mean streamlines, and one is 
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led to enquire if a simple dimensional argument will yield the result about 
the time scale. Presumably the explanation lies in the fact that the 
interval ¢ since the marked particle was released at the origin is also the 
interval of time since the turbulent eddy which surrounds the particle was 
generated at the origin x = 0. Eddies are generated with infinite energy 
and zero linear dimensions at the (virtual) origin, and then convected 
downstream, and the development of the free turbulent shear layer as 
a function of x is essentially a process of diffusion and decay of the eddies 
under the action of inertia forces. In the absence of more than one 
dimensional parameter specifying the conditions of generation of the eddies 
at the origin (such as the momentum flux in the case of a jet), the time scale 
or representative period of the eddies moving downstream must, for 
dimensional reasons, be proportional to the decay time ¢. This repre- 
sentative period of the eddies is identical with the time scale of fluctuations 
of the velocity of an element of fluid in the eddy, and so the above general 
result is recovered. 

The second interesting common feature of the results of §5 and $6 
is that 


D,(t)  [X(t)]’; (8.2) 


for both jet-type and wake-type flows the longitudinal dispersion increases 
in proportion to the thickness of the shear layer at the mean position of the 
particle. (‘The same is true of D,(t) in cases of two-dimensional flow, 
as a consequence of the proportionality between D,(t) and D,(t).) ‘This 
can be regarded as essentially a product of a dimensional argument, although 
such an argument might not be convincing by itself in view of doubt about 
whether the thickness of the shear layer at the mean position of the particle 
is really the only length available as a measure of the dispersion. (Results 
appropriate to a case in which the marked fluid particle is released at some 
finite value of x may also be obtained with analysis like that already described. 
The same forms of D,,(t), etc., as those given above are found for sufficiently 
large values of t, but the results for smaller values of ¢ depend on the position 
of release and are certainly not obtainable from dimensional arguments.) 


9. "(THE MAXIMUM MEAN CONCENTRATION OF MARKED FLUID 

As the marked particle moves downstream the statistical dispersion 
of its position increases, in the way already described, and the probability 
of finding the marked particle in unit volume located at its mean position 
decreases correspondingly. ‘This probability of finding the marked particle 
in unit volume at any given position and time, which may be termed the 
mean concentration of marked fluid, has its absolute maximum value at 
some constant value of (y/x”, z/x?) which is unknown (although for any 
flow which is symmetrical the maximum clearly lies on the centre line), 
and at a distance downstream x = X(t) which increases with ¢. It is implicit 
in the hypothesis made in §4 that the whole probability distribution of the 
displacement X(t) attains a self-preserving form, and the same will be true 
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of the displacements in the two lateral directions; thus the maximum mean 
concentration of marked fluid, C,,(¢) say, may be determined from the 
relation expressing conservation of the total amount of marked fluid. ‘The 
linear extent of the distribution of mean concentration of marked fluid in 
the v-direction is measured by D,,(¢), in any lateral direction in which the 
turbulence is of finite extent by [.X(¢)]”, and in any lateral direction in which 
the turbulence is stationary by D,(t). Hence 


C,,(t). D,(t). Lx)” = constant 


in cases in which the turbulence is of finite extent in the two lateral directions 
(as in a round jet), and 


C,,(t). D,(t). D,(t). X(t)” = constant 
in cases of two-dimensional mean flow. In all cases we have 
D,(t) « D,(t) « (X(t), 
so that the variation of C,,(t) can be written generally as 
C,,(t) «& [X(t)]}-*”. (9.1) 


If the release of marked particles at x = 0 is continuous, the distribution 
of mean concentration of marked fluid is steady and the maximum mean 
concentration (now only a maximum with respect to y and 2) is a function 
of x alone, say C;(x). For turbulence which is of finite extent in the two 
lateral directions, uniformity of the flux of marked fluid across different 
sections downstream requires 

xt = xt" for jet-type flows, 


C,,(x) oc ¢ ° 


Q? 
cl . (9.2) 
aie = for wake-type flows. 


For turbulence which is stationary in one lateral direction, the diffusion 
in this lateral unbounded direction must be taken into account, but, as 
already seen, the dispersion of a particle in this lateral direction is propor- 
tional to the thickness of the shear layer at the mean position of the particle, 
and the above functional forms of C(x) are again applicable. (‘These 
results for C/(x) can also be obtained by integrating the contributions 
from a whole set of instantaneous sources of marked fluid, released at 
different times, provided the fact that the similarity laws are applicable 
to wake-type flows only when wx and ¢ are large is employed.) 


10. TABLE OF RESULTS FOR THE VARIOUS FREE SHEAR FLOWS 
The values of the similarity powers p and q for the cases of the two- 
dimensional jet, wake and mixing layer (with zero velocity in the stream 
on one side) and the axi-symmetrical jet and wake are well known (Goldstein 
1938). Substitution in the formulae (4.4), (5.1), (5.5), (6.4), (9.1) and (9.2) 
then gives the power laws shown in table 1. 
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| at | 
| Round | Plane | Mixing | Round |} Plane | 
|} jet jet | layer | wake | wake | 
| 
Sulerian length scale L(x | | | LE 4 | 1 
Eul l I le L(x) | 5} 
Eulerian velocity scale V(x) | 1 A 0 | 5 4 
; : =. | A 
Mean displacement of particle X(t) | 4 5 f | 1 l 
Lagrangian velocity scale w(t) | 1 1 0 | 3 4 | 
= | 
; ; 
Lagrangian time scale 7(t) | l | 1 | 1 1 
Dispersion D,(t), D,(t) I 3 t | ; 4 
Maximum mean concentration | | 
| c ° ’ | 2 3 
(release of one particle) C,,,(t) = 3 2 5 | 1 3 
| Maximum mean _ concentration 
(continuous release) C,;,(x) | I 3 7 & 1 
| | 








‘Table 1. Showing powers of the independent variable (either x or t) for the quantities 
listed on the left hand side. 
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SUMMARY 

‘his paper presents the application of the methods developed 
in a previous paper (Part I, Crane & Pack 1957) to the mixing 
of two parallel streams for both laminar and turbulent flows. 
The effects of both high velocity and large temperature difference 
are treated together. ‘The method used consists in developing 
the stream function in a double series of powers of two parameters, 
the first being the Mach number and the second depending on 
the temperature difference of the streams. Analytical expressions 
are found for the terms up to the second order in the series for 
tl 


in velocity and temperature. However, when one of the streams 


(=) 


e stream function when the streams do not differ too greatly 


is at rest the analytical method is no longer sufficiently accurate, 
and for this case numerical solutions are given. 

For laminar mixing the most important effect is that of ‘ change 
of scale’, as was found in Part I for a laminar jet at large distances 
from the orifice. For turbulent half-jets the effect of ‘change of 
scale’ and the effect of the perturbation terms due to the Mach 
number of the flows are approximately equal and opposite, 
leaving the form of the velocity profile sensibly unchanged from 
that in incompressible flow. ‘This last result is confirmed by 
comparison with some experiments of Laurence (1955) on a 
two-dimensional jet at M = 0-7. Lastly, the effect of temperature 
differences is shown to be relatively unimportant even when 
these are fairly considerable. 


INTRODUCTION 

‘The mixing of two semi-infinite incompressible streams has been studied 
by Gortler (1942), by Lessen (1949), and by Lock (1951). Gortler’s method 
of analysis, applied to turbulent flow in the streams, was adapted by Pat 
(1954) to the study of laminar mixing of incompressible streams. Gortler’s 
solution is here used as a starting-point for the calculation of the mixing 
of streams of compressible fluid. ‘The functions originally calculated 
numerically by Gértler are given below in analytical form. It has been 
found necessary, to treat separately the cases (i) where the velocities of the 
streams are not very different and (ii) where one of the streams is at rest. 
This is because the convergence of the series obtained by Gortler’s method 
is not sufficiently rapid in case (ii). 


F.M I 
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EQUATIONS OF MOTION 

Let u, v be the velocity components (or mean velocity components in 
the case of turbulent flow) parallel to cartesian (x, y)-axes. Let the origin 
of coordinates be taken at the point at which the mixing begins. Let p be 
the density of the gas in the jet, and 7 the absolute temperature. ‘he 
assumptions on which the theory rests in this paper are exactly those used 
in Part I (Crane & Pack 1957) in establishing the equations of motion. 
These assumptions are that the boundary layer equations hold and that 
in the turbulent case there exists a coefficient of eddy kinematic viscosity « 
which is independent of the y coordinate. Hence the basic momentum 
equation for both laminar and turbulent half-jet mixing is equation (6) 
of Part I, namely 


— ae a ee ae Se a? Ts a (1) 
OS OCCS 0g os” Os | Cs” 
In this equation, % is the two-dimensional stream function defined by 
pu = o/dy, pu = —od~/dx; z is a variable defined by 3 = ¢ py av; 
~O 
T* = T/T, (the suffix 0 refers to a fixed state of the fluid to be defined 


later); and ¢ is a function of x. For laminar flows € = x, « = py and 
8 = n—1,n coming from the law of variation of the coefficient of viscosity pu 


with 7 (i.e. up. 2 T”). For turbulent flows €= | e(«) dx, where e(x) is 


an experimentally determined function of x defined by € = eg) e(x), €y being 
a constant; also, x =€)p) and 6 = —2. Let the two streams in their 
uniform state have speeds U,, U, parallel to the x-axis, as they cross the 
half-lines y > 0, v < 0 respectively. Let 


U,=4(U,+U,), A=(U,—U,)(U,+U,). 


The solution giving similarity of velocity profiles is expressed by 


U> on\ 2 
bb = (py Uy xl)! “g(n); a eo s 
a J 
and the equation for g(7) is 
“ STR oO" 4 Logo” — () ? 
dy | 5 SV 28s. 7s (2) 


where dashes indicate differentiations with respect to 7. 

As in Part I, ifthe Prandtl number is assumed to be unity then a particular 
integral (Crocco’s relation) satisfying the boundary layer equations of 
energy and momentum for both laminar and turbulent flow of jet type is 

lw+i1= A+Bu, 


where A, B are constants, the values of which depend upon the boundary 


conditions in the given problem, andi = e 6 i (the enthalpy per unit mass). 

Crocco’s relation has to satisfy the conditions u = U,, T = T,, p = p, at 
T ry" ry ¥ ry ’ + . 

y= + w,andu = U,,T = T3,p = pgaty = — 20. Theselead to the following 


expression for the Crocco relation: 


/ 
= 1+ : (u* —1)—w? u*?, ( 


o>) 
~— 
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where 
7° = T/T, u* = u/U,, w. = Us 2C pT: 
I+, ck ca I,-1, 
,=—_=p -= C_f,+4U- = 1,2 = ——_—_-_, 
ly =. , I, C p T; 2 i (2 i? )» h I, hi i 


The properties jg and py of the gas are taken at the temperature 7). ‘The 
quantities J, and J, are the stagnation enthalpies in the reservoirs from 
which the uniform streams may be supposed to have come. 

The equation (3) shows how the initial discontinuities of temperature 
and velocity are smoothed out when the mixing is completed. In fact, 
the stagnation enthalpy at a point in the jet where the velocity is u, obtained 
by bringing the fluid at that point to rest adiabatically, is equal to 

(1,—1,)(u— U4) 


lie a | as 


Let g(y) in (2) be expanded in the double series: 


&(y) > 


7 ) 


Y rg fy (E)wi'h'; (4) 


where € = }(y—)), 7 = b giving the locus of those points of the flow where 
the velocity is Up, and all the a,, are equal to — f except ayy which is equal 
to 2. ‘Then 


u aa . -/ 2 
sie ae S S te. 7 ark 
y% ae." 
and 
rym »] h 5 . 
(| =] 39) 5 (u* 1) Wu (>) 


‘The above expressions are valid for sufficiently small # and w>, since the 
numerical value of (w*—1)/A is less than unity. 
The boundary conditions on g’(y) are 


g(+o)=14+A, g(d)=1, g(-—a)=1-A. 


Lock* (1951) has shown for the incompressible case that any given solution 
of (2) generates an infinity of solutions, all of which satisfy the same boundary 
conditions at 7 = + «. ‘These solutions may be obtained by replacing » 
by » +c in the given solution, c being an arbitrary constant. ‘This result 
can be shown to be true for the compressible case. Each one of this infinity 
of solutions is an equally valid solution of the boundary layer problem 
because the solution obtained by replacing » by 7 +c in any given solution 
leads to a value of the y-component of the velocity on the boundary which 
differs by an amount of order U,c/,/(Re) from that of the given solution, 
Re being a large dimensionless constant. (In laminar flow Re is the usual 
Reynolds number; when the flow is turbulent Re is defined as in the 
laminar case except that the coefficient of kinematic viscosity j19/py is replaced 
by the eddy coefficient of kinematic viscosity €).) ‘l’o pick out the correct 

The author is indebted to the referee for drawing his attention to Lock’s work, 
and also for other helpful criticism. 
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solution it is necessary to take into consideration those terms of the 
Navier-Stokes equations of order 1/\/Re higher than the terms of the 
boundary layer equations. It follows from the above discussion that the 
velocity at any finite » is mathematically indeterminate from the boundary 
layer equations alone. ‘Thus 6, which by definition gives the locus of 
points in the flow whose velocity is Up, is likewise indeterminate. ‘The 
value of b can however be found from experiment. When 4 is known the 
solution to the problem is uniquely determined. 

The equations for the functions f,, are obtained by substituting (4) 
and (5) into (2) and equating the coefficients of w7"h, to zero. Thus 


Soo + 2foo foo = 0, (6) 

. : 7) 1 ; d i) ol _ 

fio t 2( foo Tio + foo S10) va {foo Foo} = 0, (7) 
dé 


foo— 1 ) 
-2( foo for Foo to) - ~ 2 A ae r = V, (8) 


ind so on. ‘The boundary conditions are 

Sook + 00) = 1+, Foo(0) = 1, Jal - 00) = a 

fh + ©) = 0, fio(9) = 0, fio — 0) = 

foul + 0) =9, foi(0) = 0, fil 00) a 
The solution of these equations will be found for two separate cases, namely, 
when A is small (say less than 0-3) and when A = 1, the latter corresponding 


to one of the streams (here the lower stream) being at rest. In this way 
many examples of practical interest are likely to be covered. 
U e 7 


Case I. SOLUTION FOR SMALL VALUES OF A 


(1) The functions f,, may be expanded in series of ascending powers 


| 
tnen 


F'(0) = 0 Fi(+2)=0, (n= 2, 3, «). 


’ n 


The functions F’, are to be derived from differential equations obtained by 
inserting the above expansion of fy) in (6) and equating to zero the coefficients 
of successive powers of A. ‘They were first evaluated numerically by Gortler. 
‘They may also be expressed in analytical form, and have been tabulated 
here from these analytical expressions (see table 1). ‘The numerical values 
show that the expansion of fy) does not converge very quickly when A is 
greater than about 0-3. For this reason the important case when A = 1 is 


considered separately below. 
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The following are the analytical forms of Gortler’s functions, in which 
4 : ’ 





2 5 ? 
O(2) =— | ed, 0,2) = Te“, 
? (¢) \7 , 2 ls) bad 
V7 [1 l ae 
—_ (3 -) 0-72521 
Fy = &*, 
F, = a - €D rw 10,4 + B, 
F, = — €0?— 300, + 8, O+ \/(2/7)O(v2E) + 3€, 
F, = 403 — (429 + 3£)@? —— — (422 + POF O — [6 ED} — (3v 3/477) O( 3€) — 


1 
4 
~ B,(E2 + )O®, — 48, EO? + (3/3/47 — 4) + (4 — B20, 4 


(F, is not evaluated since F,(0) = d, which is unknown.) 
(2) Next let 
to = GIG. 
a 
The boundary conditions are 
G10) C+ 0) — 0, = 9, 1, «1 
The equations for G,, obtained by insertion into (7) are 
GC” 426s" =P, 
G" +2£G" +26, F+2F" =0, 
andsoon. ‘The solution to the first of these equations satisfying the boundary 
conditions is found to be G) = ad), where dy, is a constant to be determined. 
In the equation for G,, put G/ = ®, v(€). Then 
: : 
y’ +2a,—4é = 0. 
Thus G, = (26 —2a,€+C)®,, 
On integration and application of the boundary 


where C is a constant. 
-1, a, = 0, and that 


conditions, it is found that C = 
ee 
G, SO, t a}. 


Phe value of a, is found when the solution for G, is calculated. In fact, 
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Gy = —207+00,(& — 3€) + OF(32 — 1) + (28, P- Bi +a, 
‘The analytical forms of the functions 


? 
Jb, +2 
where f, has the same value as before. 
G,, G,, Gy are therefore 


G, = 0, 
; , 5 oe 
1 =39,+4,, with a, = v7| = = 1-36927, 


' ree Pe 2’ ; : 
G, = -4*- (58: 4), (5) D(y 2) — FEO] — B, EY, + a, D+ 2é. 


The solution satisfying the boundary conditions on Fy is Fy = &+-c, where c is a 


On insertion of this expression for Fy in the equation for F’, it is seen 


constant. 
that c = 0 if F, is to satisfy the boundary conditions. 
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‘The expressions for G3, G,, etc., could be found in the same way, but 
the labour involved in equating the complicated analytical terms would be 
very great, for example, G, alone involves about 60 terms. When the 
above expressions for G; and G; are used, f,,, is obtained to 1°, accuracy 
for values of A up to about (0-2 if the Mach number is less than five. 


(3) ‘The same form of expansion is used for fy,, namely 
to = > A, (€)Ar", 
) 
with the boundary conditions 


H'(0) = Hi(+a)=0 (n=0,1,...). 


A procedure similar to that used above then yields the results 

Hy = ¢, = $vn(1—2/7) = 0-32204, 

H, = c,0+é&@? +4 SOO, —&, 

i Le bi EMn + FAW eet. R 2 1) qy3 
A, ny (38; a 3C}) ~S (Cy a 28,)}O®, — &(d¢ B,)®; 7 (3)® Li 
(JE — 2)D2D, — 20? & + (3/7) (032) — 
— JED} + (3 — V3/7)O — 2B, c, ED). 

(H, has not been determined.) 

As for the other functions, the expressions for the higher H-functions 
could be found exactly, but the labour would be very great. However, 
when only the functions up to Hy, are used, f), is found quite accurately 
. . / / . 
for values of A up to about 0-3. Indeed, on comparing H, +H, with the 
exact solution for f,, when A = 1, it is seen that in the range € 1:5 the 


agreement is quite good. ‘lhe maximum error in this range is about 5°... 


Cask II. SOLUTION FOR A | 


When A = 1 the lower stream is at rest. The solutions of (6) to (8) are 
required subject to the conditions 


ae Le 20) = hs fo(9) = l, fil- KL ) =U, 


with all the other conditions unchanged. 

An iterative method was used to solve equation (6). ‘he solution 
tabulated below is in agreement with that of Lock (1951) for the same 
problem. 

The next equation to solve is (7). Again a numerical method ts required. 
It is seen that an integral belonging to the complementary function is 


tie 2a Poy By the use of this integral it is possible to reduce (7) to a 


second-order differential equation. When the asymptotic form of the 
complementary function of (7) is examined it is found that (7) has solutions 


which result jin f,) remaining finite at €= + © whereas the reduced 
(second-order) equation has integrals which diverge exponentially at 
€= +0. Thus it is advantageous to solve (7) itself numerically, rather 


than this reduced equation. 
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‘Two linearly independent parts of the complementary function of (7 
with boundary conditions 
fro) =9, fio) = 9, fi,(0) = 1, 
and fo =1, fa(0)=9, fio(0) =9, 
respectively, were computed. A particular integral of the complete 
equation (7) with boundary conditions 


-/ 


ff0) = 9%, f £0) = 0, fi 0) = 0, 
was obtained. <A linear combination of the first two solutions was added 
to the particular integral to give the solution satisfying the boundary 
conditions. ‘This treatment of (7) is equally applicable to (8) and was 
used to find the values of fy,(€). Values of fy, for, fio and their first 


derivatives are given in table 2. 


‘THE EFFECT OF THE CHANGE OF SCALE 


If w> and h are both small, 
€:=E+ 5] (fow—1) d&—-4 | (foo)? 4 (9) 


where 
¢ ' {l 0 Po = 
$4>3| > \ 
arias 
When A is small, this approximation becomes 
g e £4 Ati Fy F(0)} a AF, u V?\ Fs - F(0)} T ++] ‘ 
) —wifE+2\ F, — F,(0)}+...], 
but when A = 1 the equation (9) should be used. 
‘These relations show how changes in the width of the mixing region 
due to the ‘change of scale’* depend on the difference of the stagnation 
enthalpies of the streams and on the Mach numbers of the flow. ‘To a first 


approximation these effects are independent. By considering | (fio—1) dé 
~ oO 


it can be shown that the net effect on the width of the mixing region of 
differing stagnation enthalpies is very slight when hk < 0-3 for any A. 
An example of 4 = 0-3 is afforded by two streams with stagnation 
temperatures of 300°C and room temperature respectively. ‘The effect 
of a non-zero Mach number is to decrease the width of the mixing region 


as the Mach number increases. 


CONCLUSION 
The effect of compressibility on the non-dimensional velocity profile 
may be conveniently divided into two parts. One is the effect of change 
The ‘ change of scale’ effect was defined in Part I as the difference between 
(i) the velocity profile obtained from the first term in the expansion for 4, when 
z is expressed in terms of the physical ordinate y and (11) the velocity profile in 
incompressible flow (which comes from this same term with z = 4). 
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of scale. ‘The other is the change in the velocity brought about by the 
perturbation terms f,,, fin) «+. The non-dimensional velocity profile is 
given by 


U, = foot (— SB)hfor + (— 2B) fiot ---- 


‘The perturbation terms in this equation depend on both compressibility 
and the variation of viscosity with temperature for laminar fiow. For 
laminar flows the effect of change of scale is dominant, the perturbation 
effect being negligible. As an example, for A = 1, the change in the velocity 
profile produced by putting » oc T°*® instead of » a T is at most 4%, for 
Mach numbers up to 5. Furthermore, since for A < 0-3 the effect of h 
on the change of scale is negligible, the basic cause of changes in the profile 
in compressible laminar flow is the Mach number operating through the 
change of scale. ‘The larger the Mach number the narrower is the mixing 
region. For example, when the lower stream is at rest and the upper stream 
has a speed given by M = 5, the width of the mixing region is about 
three-quarters of that for incompressible flow, the Reynolds numbers 
being the same. ‘The greater part of the contraction of the profile is in 
that part where the velocity is highest—as it is expected on physical grounds, 
since the effect of compressibility is naturally less in the more slowly moving 
parts of the mixing region. When the temperature of the upper stream is 
higher than that of the lower, that is, / (), the upper half of the profile 
tends to broaden and the lower to shrink. ‘The effect is of course reversed 
when A <0. ‘This is also to be expected on physical grounds, since 
molecules in a hotter stream have higher random velocities than those 
in a cooler one and hence tend to transfer some of the momentum of the 
faster moving parts of the jet to the surrounding parts to a greater extent 
than in the cooler stream. 

For turbulent flows the perturbation terms and the change of scale have 
etfects of the same order of magnitude. It is found on computing the 
non-dimensional velocity profile that they almost cancel each other. ‘Thus 
the form of the non-dimensional velocity profile is left sensibly unchanged 
from that obtained in incompressible flow. In order to test this last result 
against experiment a scale factor « must be used*, where € = o(y—9)/x 
if the coefficient of eddy viscosity is taken (following Gortler 1942) to be 
eyx/L. y= yo gives the locus of points in the flow at which U = U,,. 
From the equations developed above, o = (U,) L/2€))'*. ‘The theory has 
been applied to the case examined experimentally by Laurence (1955), 
namely to a turbulent jet at Mach number 0-7 entering a medium at rest. 
‘lhe observations chosen for comparison were those taken in the neighbour- 
hood of the orifice where the core of the jet was at constant velocity. ‘This 
might be expected to give the closest representation to the mixing of two 
parallel streams. ‘The results have been plotted in figure 1 with u/U, 


*'The scale factor o is not an absolute constant but may vary with the Mach 
number of the undisturbed stream, that is with w?. ‘This does not affect the state- 
ment made in the preceding sentence. 
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as ordinate, U, being the velocity of the core of the jet, and with (y)—y)/a 
as abscissa. ‘The best fit between theory and experiment was obtained 
for o = 12-7, when € = o(yy—y)/x. The fit is seen to be good except in 
the part of the profile where the velocity is lower; in this region the gradient 
of the theoretical profile is less steep than that of the experimental profile. 
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Figure 1. Comparison between Laurence’s results and the theoretical profile with 
oC 12-7, for the mixing region near the orifice when a turbulent jet of air 
issues at Mach number 0°7 into air at rest. 


Now over the whole flow the gradient of the laminar profile corresponding 
to the laminar Reynolds number of the jet is steeper than the turbulent 
profile. ‘This suggests that turbulence may not be fully developed in the 
lower velocity region, that is, the flow may be only intermittently turbulent 
there*. 

* Subsequent to the writing of this paper the author’s attention was drawn 
to the work of Johannesen (1957). Johannesen has shown that the solution /,, 
gives an accurate fit, except in the region of lower velocity, to his experimental 
results for the non-dimensional velocity profile of the mixing region near the orifice 
when a round turbulent jet issues with Mach number 1:6 into a medium at rest. 


The value of o used was 21-9. 
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SUMMARY 


Hodograph methods are applied to determine the flow at 
high subsonic and sonic velocities past two-dimensional, thin, 
symmetrical bodies. ‘The boundary value problem for the 
determination of the stream function ‘’, which in the present 
theory is a solution of ‘Tricomi’s equation, is simplified by the 
assumption of a free stream breakaway at sonic velocity from 
the shoulder of the body. A solution is obtained in terms of 
Bessel functions. 

In §§2 and 3 the flow past a wedge of small angle is discussed 
and expressions are obtained for the pressure on the nose, the 
drag coefficient and the width of the wake. A comparison with the 
corresponding results in the case of sonic velocity derived by the 
more complex analysis of Guderley & Yoshihara (1950) shows 
that the present simpler theory yields very similar values for the 
pressure over the nose. 

In § 4 the flow at sonic velocity past a profile which is a first-order 
perturbation upon a wedge profile is analysed on the basis of 
the same free streamline theory. The flow pattern is obtained 
past an arbitrarily specified body by an application of the Hankel 
inversion theorem and an expression is deduced for the drag. 


1. INTRODUCTION 





In the determination of the steady two-dimensional flow of a non-viscous 


gas past a symmetric obstacle various methods, both approximate and 
exact, are available provided the velocity of the gas throughout the flow 
field is bounded away from the local velocity of sound, either above or 
below. But in the neighbourhood of sonic velocity the approximate methods 
in general break down because of the change in type from elliptic to hyperbolic 
of the relevant differential equation at this point. 
overcoming some of the difficulties is the hodograph transformation. By 
means of this, the stream function ’ can be written as a function of the 


A powerful method of 


By 
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velocity variables and the resulting differential equation is linear. Further 
it permits a solution by separation of variables in terms of trigonometric 
or hyperbolic functions and of certain hypergeometric functions usually 
denoted by #%,,(7). Here m is a parameter, real or complex, and + = q?/q> 
where q is the velocity of the gas and q,, is the maximum velocity attainable 
by the gas when subject to adiabatic expansion. A description of the main 
properties of %,(7) has been given by Lighthill (1947). Suppose now that 
the velocity of the gas is nearly sonic everywhere and that the polar velocity 
angle @ is small. ‘Then the cartesian components of velocity U and | 


(where V = 0 at infinity upstream of the obstacle) can be written as 
U = a*—-u', 
V=v0, 


where w’ and z’ are small. Here a* is the velocity of the gas when the Mach 
number M is 1. ‘The equation for ‘’ can now be approximated by 
oY y4 aay 


=m += ua = I, 
cu” a dv’? 








where y is the adiabatic index of the gas. 


By defining dimensionless variables u = (y+ 1)u'/a*, v = (y + 1)v'/a’*, 
this becomes 
2 32y* 
+ us 


Aan: 
Cue ow 





the well-known equation of Tricomi. 

Roughly speaking there are two main methods of applying the hodograph 
transformation to the problem of flow past a body. One is to use the 
solutions (7) of the exact equations and, with analogous problems in 
incompressible flow as a guide, to write the stream function as a series 
of the form t= > a, %,(7)sinn@. In general the body shape which this 
solution gives has to be determined subsequently; the coefficients a, 
cannot be specified for a body of arbitrary shape. ‘This is a very unsatis- 
factory feature of the method but the difhculties may be overcome in the 
cases where the streamline ‘" = 0, the dividing streamline part of which 
follows the body contour, consists of portions on which either q or @ is 
constant. In general such patterns in incompressible flow are preserved 
in the generalization to allow for compressibility. ‘Thus, for example, the 
problem of flow past a wedge when the streamline breaks away from the 
shoulder with constant velocity equal to the velocity at infinity upstream 
can be solved exactly for any wedge angle and any velocity at infinity, 
provided M, <1. The suffix 1 will in future refer to conditions at infinity 


upstream. 

An alternative method, much favoured by American workers in this 
field, is to use the 'l'’ricomi approximation and to set up the problem directly 
in the hodograph plane. An essential requirement in this approach is the 
determination of the singularity in the hodograph plane corresponding to 
the velocity at infinity. A singularity is clearly located here, since all the 
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streamlines both originate from and return to this point. ‘The singularity 


is consequently ‘of doublet type’ and this may also be surmised by con- 
sideration of the corresponding singularity in incompressible theory. 
However, while this singularity has been successfully identified by a number 
of authors, the same sort of difficulties arise with the other boundary 
conditions as previously. For a given boundary in the physical plane 
it is not in general possible to determine the relation holding between 
the velocity variables on this boundary and this means that the shape of 
the streamline ‘= 0 in the hodograph plane cannot be predetermined. 
However, as before, when either the magnitude or direction of the velocity 
vector is constant on ‘Y= 0, a solution of the problem is forthcoming. 

Sections 2 and 3 of the present paper are concerned with discussion 
of flow past a wedge of small angle when the velocity at infinity is subsonic 
or sonic. For such a problem in incompressible flow when the velocity at 
infinity is g,, the standard methods of Kirchhoff and Helmholtz lead to a 
solution in which the flow breaks away from the shoulder at g = q, and the 
free streamline which starts at the shoulder retains this constant velocity 
to infinity downstream. ‘Throughout the flow field g <q,. It is not 
desirable to obtain a similar solution for the case of compressible flow 
when the velocity at infinity upstream is subsonic because experimental 
evidence shows that, for reasonably high values of ,, the gas is accelerated 
up the wedge side until it attains sonic velocity at the shoulder. We shall 
therefore select a flow pattern which exhibits this feature and the appropriate 
solution of ‘Tricomi’s equation will then be found. ‘The following model 
will be adopted. ‘The gas comprising the dividing streamline Y= 0 
starts from infinity upstream and moves in a straight line up to the stagnation 
point at the tip of the wedge. It then accelerates up the wedge side, reaching 
sonic velocity at the shoulder. ‘his velocity is retained until the streamline 
becomes parallel again to the free stream. It then remains straight, 
decelerating from sonic velocity to the velocity at infinity. We note that 
the streamline Y = 0 comes within the general category we have mentioned 
in that it consists of portions along which either the magnitude or the 
direction of the velocity vector remains constant. 

‘This model has the advantage of yielding the relatively simple analytical 
solution developed in the next section while retaining the essential physical 
characteristics of flow past a wedge. In practice it might be expected that 
the Mach number of the streamline separating from the shoulder will be 
somewhat greater than 1, but the effect of this, particularly in the region 
upstream of the shoulder, will be small and may be neglected. A further 
advantage is that this is the type of model which Roshko (1954) adopted 
in studies of flow of incompressible fluid past bluff bodies. ‘The fact that 
Roshko obtained such good agreement with experimental results using 
this ‘notched hodograph’ model suggests that this is a good basis on which 
to attempt a solution. Of course, in the limit as M,— 1, we obtain the 
problem of the wedge in a sonic stream with a free streamline at sonic 
velocity extending from the shoulder to infinity. As has been mentioned, 
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this particular problem can be solved exactly, that is, in terms of the exact 
differential equation for the stream function as distinct from the transonic 
approximation, but the leading term of this exact solution for small wedge 
angles is the limiting case of the problem considered here as M, — 1. 
The expression for the stream function coincides with that of Imai (1952) 
who used a different approach. Imai based his solution on the incom- 
pressible flow and then replaced powers of the velocity by associated Bessel 
functions. His results, and those obtained by Mackie & Pack (1955) are 
discussed in § 3. 

This limiting case is important because it gives a fairly simple represen- 
tation of a wedge in a sonic stream and therefore it can be compared with the 
solution of Guderley & Yoshihara (1950). ‘This solution was derived to 
represent the flow of a gas for which M, = 1 past a wedge profile of small 
angle up to the limiting Mach wave emanating from the shoulder. Down- 
stream from the limiting characteristic any disturbance introduced will not 
affect the flow before it and consequently the continuation of the solution 
can be performed by the method of characteristics in a way which is 
determined by the subsequent profile of the object (parallel sides, diamond 
shape, etc.). When M, = 1 anew difficulty is encountered in the determina- 
tion of the correct singularity in the hodograph plane. Whereas a subsonic 
singularity is isolated, occurring as it does in the elliptic region of the 
hodograph plane, the singularity is now at the origin and will be propagated 
along the characteristic in the hyperbolic part of the plane which starts 
from the origin. ‘his characteristic maps into the limiting Mach wave in 
in the physical plane and care must be taken to choose a singularity which 
dpes not map this characteristic in the hodograph plane to infinity in the 
physical plane. ‘The correct determination of the singularity and the 
subsequent addition of non-singular terms to satisfy the boundary conditions 
is a matter of considerable complexity and to keep the analysis manageable 
Guderley & Yoshihara found a certain amount of asymptotic approximation 
unavoidable. ‘Thus it is of interest to compare the mathematically much 
simpler model of the ‘free streamline’ theory with the solution of Guderley 
& Yoshihara. The pressure (and hence velocity) distribution on the wedge 
side shows very little difference between the two solutions. A specific 
comparison is made in § 3. 

The fact that the discrepancy is small is not surprising on physical 
grounds. For both solutions describe a flow with a stagnation point at 
the tip and sonic velocity at the shoulder of the wedge, while the flow at 
a great distance from the wedge is uniform with M=1. ‘The difference 
occurs in the form of solution immediately downstream of the shoulder 
and the upstream influence of this is of necessity small in a near sonic stream. 
Thus the principal difficulty in the work of Guderley & Yoshihara, which is 
the correct determination of the flow between the sonic line and the limiting 
characteristic, does not arise here since the whole physical plane (exclusive, 
of course, of the wake) is mapped into the elliptic part of the hodograph 
plane bounded by the sonic line. ‘The absence of a supersonic region 
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means further that no limit lines can appear. ‘The occurrence of limit 
lines is an ever present danger in the hodograph theory of supersonic flow 
as their images are not readily detected in the hodograph plane but they 
render the transition to the physical plane meaningless. 

The satisfactory representation thus obtained of the flow pattern 
upstream of the shoulder justifies an extension which is carried out in § 4. 
A profile is considered which is a first-order perturbation on a wedge 
profile. ‘he corresponding flow past this profile can be obtained by an 
application of the Hankel inversion theorem. An unusual advantage of 
this solution is that it can be obtained for any given first-order perturbation 
in the physical plane and does not involve the a posteriori determination 
of the bounding streamline. An analytic expression for the drag on such 
a body is given in the general case. 


2. ‘THE WEDGE IN A SUBSONIC STREAM 
Before proceeding we shall summarize briefly the relationships holding 
between the physical variables. It should be remembered that in many cases 
these relations hold only within the limits of the transonic approximation. 
We have already defined u and v. ‘The pressure and density are p and p 
respectively. ‘The Mach number ™ is given by 


1— M?= u. 


The polar angle 6 is simply related to v by means of the equation v = (y + 1)é. 
The equation of continuity and the condition for irrotational how may be 
written respectively as 


uu,—v, = 0, u,+v, = 0. 


Inverting these, we get the hodograph equations 


uy, —x, = 0, ty, +7, = 0. (1) 


‘ 
These lead immediately to 'T'ricomi’s equation 


\ 


) uu 


+ UV wy = 0, 


in which u > 0 represents subsonic flow. ‘Thus y satisfies the same 
differential equation as ‘Y’ and it may be shown that ‘lV’ is a simple multiple 
of y. Similarly ®, the velocity potential, is a multiple of v and in much of 
what follows we shall work with x and y as dependent variables instead 
of ® and V’. In subsequent work with the flow nowhere supersonic we 
shall find it convenient to use a new variable r defined by 


r= dy, (2) 
Tricomi’s equation may then be solved according to the usual procedure 
of separation of variables and simple solutions are obtained of the type 


y = rider“ | (Ar), 


¥ 
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where @_ ,,(Ar) is any linear combination of Bessel functions of order } 


and A may be real or imaginary. 

The flow pattern described in the previous section is shown in figure 1. 
The figure is largely self-explanatory. Because of symmetry only the flow 
in the upper half-plane y > 0 need be considered. The dividing stream- 
line Y = 0 is the line EOBCD. Along EO the velocity decreases from its 
(subsonic) value at infinity to zero at O, the tip of the wedge. Along the 
wedge side OB we have v = vy» = (y +1)8 where 4 is the semi-angle of the 
wedge. On BC the flow is sonic and along CD, v = 0, while the velocity 
decreases from sonic to its value at infinity upstream. 


ere, ee ee 








5 Figure 1. Physical plane. 


The boundary value problem is now set up in the hodograph plane 
which is shown in figure 2. OBCDEO is the line y = 0 and two other lines 
of constant y are also sketched. ‘There is a singularity at the point where 
D and E coincide and in addition we must have 


y=0 onv=0, uw2J, 
v=0 onv=v%, u2>d, 


y=) tone =U, Wee Uo» 


The stagnation condition at the tip gives, in accordance with the linearization 
principle, 

y=0, x=0 asu-> o. 
Finally if the wedge is of unit length we must have 

¥=1 atv= UV, Uz Q. 


It should be noted that the problem now formulated is very similar 
to that considered by Cole (1951) in which the boundary condition along 
u=0 is not y=0 but dy/ou=0. This results from the condition he 
imposes that the sonic line be straight and normal to the direction of the 
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3 flow at infinity. A solution of the present problem may be obtained by 
a procedure similar to that carried out in §4 of Cole’s paper. It turns out 
that the required solution is of the type [’,, in Cole’s notation which he 





V : Cy 3) a 
rejects. ‘This is 
ag ° 
a (| sinhA(v,— v) 
Ss , = Aris | ——__— J 2 Ar J D Ar A dn. 3) 
. Ba },  sinhAv, Us(Ar)J y)3(Ary)A di (3) 
e Here A is a constant which will be fixed by the condition that the wedge 
is of unit length, ry, is also a constant and is the value of r corresponding to 
u = u,, the value of wu at infinity upstream. 
Vv 
4 
V=Vo 
B- : ——oO 
y= 
y=O , = ae 
ff \ 
Yi } 
( \ / 
Me! % ; 4 
Se — 
C D,E y=O 
Figure 2. Hodograph plane. 
Equations (1) now enable the corresponding x coordinate to be found. 
After a little analysis this is seen to be 
** cosh A(vy — v) 





a A(3)! 372/3 | 


I _a/3(Ar)J y)3(Ary)A aA (4) 


sinh Avy 


The arbitrary constant which appears in this determination of wx is easily 





shown to be zero because of the stagnation condition v = 0 as r-—> %. 
Since x = 1 at the shoulder where v = zw, and r = 0 we obtain from (4) 
A(3)'822/3 -> AUST, (Ar) 
l= ———- | Or 7. 
r(4) Jy sinhazry 


On expanding the Bessel function, setting Av» = ¢ and interchanging 
orders of summation and integration we get 
Ao" < (—Tpe* | See ad aie ; 
- TG) ,&, 2 813 IT 4) |, sinhe 


vo ) ) o 





1 


Making use of a well-known representation of the Riemann zeta function 
C(s), we obtain finally after some algebra 
3 x (—1)"r2"T'(n + 2)(22" +93 — 1)f(2n-+-2) 


] 
1-—4 ¥ Dit ll! 15 
7 "a (3) Pay; 22"n l u“" sl ( ) 








This equation determines the constant A. 
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In the subsequent work a series representation for x on the side of the 
wedge will be useful. ‘lo obtain this we express cosech Avy in its partial 
fraction expansion 

sesthvtes «ae ott 
cosech Az 0 a e < Pag SET ayy PERT 
Avy ni (Avp)* + (n7)° 
Substituting this value in (4) with v = vp» we finally derive, after some 
manipulation, the following series 


ioe 


: ! . ~ y ( 1)"nl 2/3 (=) K, 3 (=) ( , (6) 
©0 ( Co n= 1 a) a) ) 


for 0 <r <n, and 


ills ee y (-—1)"nlyg (=) K_ sis (=), 
l 


aT) % 


0 n 





-— 
“I 
— 


for 0. <7, <¥. 

A further series representation for x can be obtained for the portion of 
the streamline parallel to the uniform upstream flow. In this case v = 0 
and we have to use the partial fraction expansion of cothAvy. We obtain 


for OD or <r 
To summarize, we have the integral representations (3) and (4) of the 
coordinates together with representations (6), (7) and (8) for x as a function 
“of r on the straight portions of the bounding streamline. In each case the 
constant A is given by (5). 
We are now able to calculate the drag coefficient C,. ‘This we define 
in the usual way as 
’ ae ] T? 
Cy = Did, U;, 


where D is the drag on the upper surface of the wedge. ‘The local pressure 
coefiicient is defined by 
. = 1 T2 ( 
C, = (Pp —Pp,)/3p1 UT. (9) 
Then for small 6 we have 
om 
Cy =o "és ax, (10 
J 
the integral being evaluated along the side of the wedge. According to the 


linearized theory it can be shown that 


Hence 
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‘The procedure is now to substitute the series representations (6) and (7) 
for x. ‘he expression for C,, is obtained as an infinite series. Details of 





the algebra are suppressed but the final form is 
» +1)'18C cg ‘“ 
(v4 D) ‘Cp = 23/3: I ‘\e = 1)" 
>” 3 8 Ce eid 7) ! 
(7 + 8)(22" +53 — 1)C(2n + 2)(37r,/2vy)?" > —2(3r,/2v_)?%. (12) 


‘This may be written in terms of M, through the relation 
r, (1-—M?? 


3 


0 


For values of the free stream Mach number little different from unity 





the drag coeflicient may be expanded in powers of 1— Mj. ‘The leading 
terms of such an expansion are 

C 1-898°° 3 =25(1— M2) 0-491 — Me) (1 — M;)* 13 

D (y + 1)!3 ~ yt) e S13(y + 1)73— y | §7 3(y + 1/133 ? (13 


‘This expression may be compared with that obtained by Cole. As is 
to be expected, there is a fixed non-zero limit of Cp as M,-> 1. This value 
is discussed in more detail in the next section. 

In figure 1 the point C is the location of the end of the sonic line and 
the point where the wake becomes parallel to the free stream. ‘The 
x-coordinate of C is obtained by letting + — 0 in (8) and has the value 

ve Ag) ret e pa 5 = y" Pa i 
aT) XY < ) ‘ | « 0 \ al) / 

Finally we can calculate H, the semi-width of the wake. ‘This is given 
by H = 8+h where we have to evaluate 


h=(y+1)y" L dx, 


the integral being taken along the sonic line BC. Alternatively we can 
find H by taking a control surface consisting of the streamline EOBCD 
and a second streamline at a large distance from the wedge. ‘The momentum 
principle then gives 


Pi H=6 | p dx t pH —9). 
The first term on the right-hand side is the force in the x-direction 
exerted by the wedge OB and the second term that exerted by the streamline 


BC on the fluid. p, is the pressure corresponding to sonic velocity. By 
means of (9), (10) and (11) we obtain 


H = 8+My+1)Cp/u;. 


This shows that the width of the wake is infinite when the Mach number 
of the free stream is lL. 
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3. [HE LIMITING CASE OF THE SONIC FREE STREAM 


We now consider in more detail the limiting case of the solution discussed 
in the previous section when the Mach number at infinity tends to 1. The 
point C of figure 1 now moves to infinity and we have the flow of a gas with 
M, = 1 past a wedge where streamlines at constant (sonic) velocity break 
away from the shoulder and extend to infinity downstream. ‘This is the 
compressible analogue of the Kirchhotf—-Helmholtz flow past a wedge in 
the theory of incompressible flow. Although we can solve this problem 
for compressible flow in terms of solutions of the exact hodograph equations 
the relative simplicity arising from the use of Bessel functions means that 
information can be obtained with very much less labour although at the 
cost of whatever loss of accuracy is inherent in the use of the transonic 
approximation. In this connection it should be noted in particular that the 
approximation enforced by the stagnation condition at the tip of the wedge 
is an unsatisfactory feature. 

Certain of the results we shall derive may be obtained by letting 7, > 0 
in formulae of the previous section. It is necessary, however, to exercise 
some caution in performing this operation because of the special nature 
of the singularity in the hodograph plane corresponding to a free stream 
Mach number M, = 1 to which reference has already been made. For 
example, if we let 7, ~ 0 in (3), we obtain 


-© sinh A(vy—v).A 


4 sinh Av, 





4/3 
y = Br | J yg(Ar) dd, 

where Bisaconstant. Itis easy to see that the necessary boundary conditions 
on r = 0 and w = wz are satisfied by this expression but it is not clear that 
y = 0 when v = 0 as the integral does not converge here. It is possible to 
rewrite this integral as a contour integral and then as a series which gives 
an analytic continuation of the solution for v = 0. However, it is more 
natural to start from the series form as this may be written down immediately 
from the analogous incompressible problem, and then to express this as 
a contour integral. ‘The contour integral thus provides the link between 
the series and real integral forms of solution. This approach also obviates 
certain convergence difficulties associated with the derivation of the 
expression for the drag. In the series form, the solution of this free stream- 
line model has been discussed by Imai (1952). It will be seen that an 
approximation made in Imai’s computation is in fact unnecessary. 

If we make use of the standard methods of Kirchhoff and Helmholtz, 
the usual procedure in the hodograph plane in incompressible flow leads 
to the result 

Y =k" ¥ n(q/q,)'"" sin(n7O/8). (14) 
n=1 

This solution represents flow which is uniform at infinity with velocity q, 
and which has as a zero streamline a wedge of semi-angle 6 whose axis is 
parallel to the direction of flow. Downstream from the shoulder of the 
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wedge the streamline is one of constant velocity g,._k” is a (positive) constant 


1 scale factor which is generally chosen to make the length of the wedge 
e unity. 

h As previously observed, equation (14) can be generalized to include the 
k | solutions of the full hodograph equation and the physical problem just 
e described can be solved completely for any subsonic or sonic velocity q, 
n (Mackie & Pack 1955). However we shall restrict attention here to 
n performing the same operations with reference to the solutions of Tricomi’s 
s equation. ‘This consists of replacing g’*”” where it occurs by 7!’?Ky)(nz7r/vp), 
t where vj = (y+1)6. Accordingly we can write 
Y — . r'3 K,/3(n7r/v9) sin 








b Do hat AIST ways a 
<, 3 Ki 3(nzr, |v) Vo 
Now if r; = 0 a considerable simplification results, corresponding to 
the case when M,=1. Since lim z'*K,,(z) is a non-zero constant the 
. expression becomes 


y=R' > n*r3K, (nzr/v9)sin(n7v/ v9), (15) 
n=! 
where again the precise value of k’ is as yet unimportant as it will be fixed 
later. 
Cs-4¢ 
> ee 
PA 








Figure 3. ‘The complex p-plane. 


This series converges for r > 0, that is for M 1, but not when r = 0. 
To examine (15) in more detail we rewrite it as 


y=JI(W) = —4(R’) > n88r'8Ky3(nar/ve'". (16) 

We now express W as a contour integral in the form 
ea ees =e 
W = a> | v4/87/3K/4(var/v9)sin! vare’"™—-" 9) dy, (17) 


aw JC 

C is the contour in the complex v-plane which goes from —i to 70, 
indented at the origin as shown in figure 3. ‘The equivalence of (16) and (17) 
can be verified by completing the right-hand semi-circle and using the 
asymptotic properties of Ky3(z) combined with the Cauchy theory of 
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residues. ‘lhe integrand has a branch point at the origin but is one-valued 
in the whole plane cut along the negative real axis. We shall also introduce 
at this stage the contour C” (figure 3) at every point of which &(v) > 0, 
but which is such that it lies to the left of all the singularities of the integrand 
tor which A(v) > 0. 





If now we write W = W,+W, where W, comprises the part of the 
integral from () to 700, we can rewrite (17) as a real integral. In fact we 
obtain 

4 Rp? f p217/344/3,1/3 Ke in/2 —tr(1—v/%,) 

W = y Lent t4/87 8K g(e“tarr [vg )e o) — 


—e 34873 K) (e~ tarr/v))e"-"' sinh tz dt. 
Making use of the formula 
a 7/26) — 2-1/2_f6,—t7/6 - pict /6 = 
Kyjg(e" 3) = 378 Parte I _49(3) — 'S 13(8) §; 


we obtain after some algebra 
y= JI(W) = hk'n | 733743 sinh ta(1 —wv/vp)sinh |! tad 1/,(tar/v) dt. 


/ 0 


With the substitution ¢ = v,A/z, this becomes 


y = §R'x(v/7)378 | A483 sinh A(vp — v) sinh! Avg J 1/(Ar) dA. 
Comparison with (3) shows that this is precisely the same expression 
as would be obtained by letting 7, > 0 in the work of the previous section. 
We can now write 





’ yo | A*3 sinh A(vy — @) sinh™! Avo J 4)3(Ar) dA, (18) 
wher« | 
a2 33 
k = (3)"8 asa - 19) | 
TON we 


‘I'he value of k is deduced from the limiting cases of (3) and (5) as 7; > 0 
but may also be obtained independently. 

Equations (18) and (19) together give the solution past a wedge of unit 
length with a free streamline breaking away from the shoulder at sonic 
velocity. 

Since this result forms the basis of the work in the following section 
we shall discuss the solution in more detail. First we note that (15), the 
imaginary part of (17), and (18), are all representations of the same solution, 
subject to the suitable correlation of the constants k and k’. However (15) 
does not converge when r = 0 for any v while (18) does not converge when 
v =( for any yr. By a suitable choice of C’, the contour in the complex 
v-plane, we can make the integral expression in (17) convergent everywhere 
except at r = 0, v = 0 where it must have a singularity. One way to do 
this is to choose C’ to be the contour consisting of the pair of straight lines 
arg(v) = +(s7—e) where 0 <« ‘a, indented if required at 0. ‘The 
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distortion of the contour from C to C’ is easily shown to be valid for some 
range r > 0, v > O and the extension to other values of r and wv follows 
from an appeal to the principle of analytic continuation. 

The value of the x-coordinate in terms of r and v can be obtained by the 


methods of the previous section in either of the forms 


¢ = —(3/2)'3k'r?3 S n*8K,.(nar/vy)cos(nzv/vp), (20) 


x= (3 2) 8k | A4’3 cosh A(v, — v)sinh—! Av, J »3(Ar) dA. (21) 
The pressure can be obtained similarly. However, since these expressions 
again show non-convergence either when r = 0 or ¢ 1), we shall indicate 
how a rigorous formulation of the drag coefficient can be obtained by means 
of the contour integral representation. 
We have 
L / 
Cy, = 26y+1)' | udx dx, 


If now we replace y, by the appropriate contour integral derived from 
(17), we can invert the order of integration and obtain expressions such as 
[243K (2)|0 , [s?3Ag,(s)]o . ‘The reason for replacing C by C’ now 
becomes clear. For since we have performed the integrations with respect 
to u (or r) first, we must have z in the above expressions with a positive real 


part and this can only be done if A(v) Q on C’. After some algebra we 
obtain 
Saver fv, \?3 di y2/8 dy) —1 
Cp = 2 FGI te) hh sees eee 
\ 7 cr SIN 17 cq SIN viz 


C’ can now be replaced by C and the resulting contour integrals are 
standard forms for functions associated with Riemann zeta functions. 
Further algebra now reduces the expression for C,, to exactly that obtained 
when 7, is put equal to zero in (12). Numerically we have 


Cy, = 1:89893(y + 1)-18, (22) 


Some explanation should be given regarding the discrepancy between 
the factor 1-89 appearing in (22) and the numerical results obtained by Imai 
(1952) and by Mackie & Pack (1955). ‘There is an error in computation 
in Imai’s work, and in addition he replaces K,,{z) by its asymptotic form 
in assessing the drag: whereas in fact this is not necessary as the expressions 
can be integrated exactly. Mackie & Pack obtain the drag coefficient for 


small 6 as the leading term of the expression for general 6 when the velocity 
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‘igure 4. Chordwise pressure distribution. 
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is sonic. Although the details are not given, the analytic expression is the 
same as that given by the leading term in (12) as 7, > 0 but the numerical 
factor is subsequently computed incorrectly and this accounts for the 
discrepancy. 

Comparison with the more elaborate computations of Guderley & 
Yoshihara (1950) based on a solution valid up to the limiting Mach wave 
shows the pressure obtained in the present solution to be slightly higher. 
Figure 4 shows the variation of pressure along the wedge face compared 
with that calculated by Guderley & Yoshihara and the variation is seen 
to be slight. For application in the following section it is desirable to have 
some estimate of the value of r along the wedge face and this is shown by 
plotting r against x in figure 5. 
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4. ‘THE PERTURBATION PROBLEM AND ITS SOLUTION 

We now consider the flow past a symmetric wedge-type profile of 
arbitrary shape provided this is such as may be represented in the physical 
plane by a small perturbation upon the wedge with straight sides. As in 
the previous section, the Mach number of the flow at infinity upstream is 1. 

‘The pattern in the physical plane is much as in figure 1 except that the 
line OB is replaced by a curved line lying close to it. We shall take the 
semi-angle at the tip of the new body to be 6 as before but the slope of the 
profile at the shoulder will now be 8’. The equation of the boundary is 
taken to be vy = 6x+eF(x). In the subsequent work we shall neglect terms 





Figure 5. Chordwise variation of r. 
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which are O(c”). The streamline Y = 0 downstream from the body will 
again be presumed to be one of constant velocity which is the velocity of 
sound. 

We now set up the boundary value problem in the hodograph plane as 
shown in figure 6. (u,v) will have a given singularity at the origin DE 
and vanishes on EO and on BD. B is the point (0, v,) where uv, = (y + 1)d’. 
Further, y = 0 on OB where OB is some unspecified curve into which the 
body profile is mapped. Its equation is of the form v = v)+ef(r), where 
we recall that r is related to u by (2). Finally, there is the stagnation 
condition y = 0, x = Vasr> om. 


Figure 6. The hodograph plane for the perturbation problem 


/(v) can now be found. Since OB is a streamline it follows that along it 





a ee ee. ee f(r) 
> o+ef (x) mol} gi 4 = o+ y4 1° 
Thus 
[F'(x)], = g(r,v) = (y 7 1)" f(r) On VU = UT ef(r). (23) 
We now write 
x= x* +x, V=MF+EY,y, 


where x*(r,v), y*(r,v) are the values of x and y in the unperturbed flow 
obtained in the previous section for given values of r and v. If this is 
substituted in (23) we obtain 


f(r) = (y +1) 'fx*(r, 2») } + O(6). (24) 


The function x*(r, v9) is given by (20) or (21) and is shown graphically 
in figure 5. ‘The shape of the streamline in the hodograph plane is now 
given to the first order in «. ‘This result is important for it enables the 
problem to be set up in the hodograph plane directly from a given con- 
figuration in-the physical plane. ‘Thus it avoids the disadvantage of the 
a posteriori determination of the boundary streamline to which reference 
was made in §1. 
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‘The boundary condition on OB requires that y(r, v) = Vonv = vy + f(r) 
‘That is 
y*{r, Up tef(r)} +ey,{r, Uy tef(r)} = 0. 
Since y*(r, vy) = 0 we must have 
yA, %) = —f(r)[ey*/ov], — »,- (25) 

‘his function on the right-hand side, which is known when the shape 
of the profile is known, we denote by H(r). 

In seeking a solution of the perturbation problem we must therefore 
find a function y,, which vanishes on v = 0 for all u, onu = Ofor0 < v < v%, 
and which takes the value H(r) on v = vp. ¥,, is non-singular at the origin 
as the doublet-type singularity necessary for the flow is already contained 
in y*, Finally y,, must not upset the stagnation condition at x = 0, v = 0. 
All the conditions except (25) are immediately satisfied by a solution of 


the form 


~ft 


y= | G(A)r#8J,.(Ar)sinh Av da. (26) 
The condition (25) is satisfied if 
| G(A)sinh Avy 137, (Ar) dA = H(r), 
and it follows from the Hankel inversion theorem that 
A , 
- f 23) me : D7} 
ey ae | H(t)t?J,).(At) dt (27) 


0 if) 


G(A) 


‘The complete solution of the perturbation problem is now y = y* +ey,, 
where y* is given by (18) and vy, by (26) and (27). ‘The corresponding 
value of x follows from equations (1). It is x = x* + ex,, where x* is given 
by (21) and 
x, = —(5)"* | G(A)r?3J _9/,(Ar)cosh Av da. 


yp 


“i 


An expression is now found for the drag coefficient on the upper halt 


of the profile. In estimating this, allowance must be made for the variation 
in slope of the profile and for the fact that the total length is no longer 
necessarily unity. ‘lhe natural extension of the definition is 


7 ) 


0 res 
Cy = | C,,[8 +(y + 1) f(r)] ds / ds, 
’ fos 
where the integration is taken along the side of the profile and the limits 
indicate that it is carried out from the tip to the shoulder. Now on the 
profile 
Y= x" + ex, y =€Yp- 
Hence ds = dx}1 + O(e?)! and consequently we can replace ds in the integrals 
by dx and write 





=() . { ox ox, | 
| C,f8+ey+1) f(r} § a +e a4 > dr 
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‘The terms independent of €« in the numerator and denominator are 
respectively the drag coefficient of the original straight wedge given by (22) 
which we shall here call Cj, and its length which is unity. ‘Then if 
Cy = Ci, +€Cp, we have, after some simplification, 


. 36 CY, \ 
Cp) =e (y 4 De! rH(r) dr ere, | (Ze) 7 dr 4 


r(3"7Ch| r » (Ze) dr. (28) 


Inspection of (28) shows that three effects contribute to the change in 
the drag coefficient. ‘The first term represents the effect due to the variation 
in slope of the profile in that the angle between the direction of the free 
stream and the direction in which the pressure acts (namely, normal to 
the surface) varies along the body. ‘lhe second term is an estimate of the 
effect of the change in the magnitude of the pressure along the body due to 
the perturbation of the wedge profile. ‘lhe third term is a correction due 
to the change in length of the body and is necessary since C),, is defined 
as a dimensionless quantity. We now recall that H(r) is by definition 
y,(7, Up). It follows that the integrands of the three terms in (28) are all 
of the same order in 6 and hence that the terms themselves are in descending 
order of magnitude for small 6 since CZ, is proportional to 6°%. It is of 
interest to note that an assessment of the most important term can be made 
without obtaining the explicit solution. For from (24) and (25), A(r) 
depends only on the given perturbation in the physical plane and on the 
solution of the unperturbed problem. ‘Thus for a given profile, the leading 
term in the drag coefhcient, C,,,, can be quickly evaluated by a numerical 
integration, that is, by evaluating the integral | C}(dy/dx) dx* along the 
wedge, where C} is the pressure coefficient for the unperturbed wedge 
profile and dy/dx is the true slope of the new profile. 
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REVIEWS 


Rheology : Theory and Applications, Volume I, edited by F. R. 
ErricH. New York: Academic Press, 1956. 761 pp. $20.00. 


The volume under review is the first of three which will attempt to set 
out the current position of rheological theory and practice. ‘The scope of 
the individual volumes has been given by the editor as follows: ‘.. .the 
reader will find two introductory chapters, one from the physicochemical 
and the other from the physics and engineering angle, followed by five 
chapters on various phases of the deformations of solids. The paper on 
flow under high pressures leads to those on the mechanism of liquid flow, 
large elastic deformations, viscoelasticity, and melt flow. Four chapters 
on the basis of the rheology of disperse systems and one on acoustic 
responses of liquids complete this part. 

Volume II will open with an integrated survey which will serve to link 
the fifteen chapters, woven through the three volumes of the book, that 
deal with various fields and aspects of linear viscoelasticity. Volume I] 
will continue with relaxation theory and three chapters on experimental 
techniques; then there follows the series of chapters on special types of 
materials or behaviour such as the relaxation of polymers, the rheology of 
elastomers, glasses, cellulose derivatives, and fibers; it will include also 
chapters on concrete and on seismic measurements. 

Volume III will contain more specialized chapters, on crystalline and 
on cross-linked plastics, polyelectrolytes, latexes, inks, pastes, and clay. 
This part will conclude with a series of technological articles on lubri- 
cation, spinning, molding, extrusion, and adhesion and a survey of the 
general features of industrial rheology.” 

No artificial attempt has been made to force a unity on the subject 
matter for which it is not yet ready, and the various authors of the different 
chapters have had considerable freedom in the manner of presentation. 

Although the behaviour of materials under stress has been one of the 
most important characteristics of the human environment, from the 
earliest prehistory of man, it is only in the last thirty years that ‘ rheology ’ 
has grown up as a separate field of study. Much of the subject matter of 
rheology originally fell within the control of the craftsman, who developed 
considerable skill in the subjective assessment of the state of the materials 
with which he worked. Certain aspects, such as the behaviour of elastic 
solids when undergoing small deformations, and the flow of Newtonian 
liquids, were developed early as branches of physics and applied mathe- 
matics. ‘This was possible, because they combined practical importance, 
in a wide range of applications, with a formal simplicity which rendered 
them amenable to mathematical analysis. 

The last thirty years has seen a similar trend of historical development 
in most researches on substances whose mechanical behaviour has practical 
and theoretical interest. The first step is the replacement of the sub- 
jective judgment of the craftsmen by a simple test designed to give 
quantitative results, although these are usually only of the pointer reading 
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type. Secondly, attempts at precise correlation of these tests with perfor- 
mance can lead to modified test methods better able to predict the behaviour 
of the material in use. Up to this stage it can barely be said that a scientific 
analysis of the problem has been commenced. At the next stage, however, 
the problem becomes a physical study in which the relations between 
stress, strain, rate of strain, and the dependence of these relations on the 
physical variables, may be undertaken. Mathematical analyses are 
required, to express the relations in the simplest manner, and to permit the 
solution of particular problems. Finally, the macroscopic properties are 
regarded as the outcome of the structure of the material, down to atomic 
and molecular dimensions, and theories can be constructed to give quanti- 
tatively the macro-behaviour in terms of molecular and atomic magnitudes. 

This course of development is almost invariably followed in research on 
traditional materials. With new materials made through the intervention 
of science, the early stages can be by-passed to some degree and the 
investigation may, from its inception, be studied on a strictly scientific basis. 

Different materials have now reached widely varied stages of study, 
determined in part by the complexity of the problems involved, and in part 
by the research effort available. ‘This is well illustrated by the contrast 
between metals and polymers on the one hand, and soil, concrete and 
related building materials on the other. Although the rheological 
properties of the latter are of the greatest practical importance, research has 
gone little beyond semi-empirical measurements. Much of the work on 
metals and polymers is now concerned with the relation between properties 
and micro-structure, and real progress has been made in developing 
adequate theories. 

The uneven state of the subject is already reflected in this first volume. 
It will become even more noticeable when the second and third volumes 
have been published. Six of the seventeen chapters in Volume I are 
concerned almost solely with polymers either directly or by implication. 
This is not an unfair proportion in relation to the amount of published 
work—on the contrary it represents in some instances a difficult exercise 
in selection and compression (e.g. in Frisch and Simha’s chapter on “The 
viscosity of colloidal suspensions and macromolecular suspensions ’’). On 
the other hand it unavoidably fails to convey the relative practical 
importance of rheology in other fields. 

It remains to consider the purpose the whole work will serve and the 
extent to which the first volume fulfills it. The preface states that “‘ The 
contributions ... are so planned that scientific workers are introduced 
to well-demarcated areas of rheology through introductory and des- 
criptive material which then leads into integrated surveys of the present 
knowledge in these areas”’. ‘The readers whom the book will serve primarily 
are the experimental physicists, chemists, and to a lesser degree biologists 
and others, who require to make rheological measurements or use the 
results of such work. ‘The small group of applied mathematicians who 
have been responsible for much of the modern mathematical theory of 
rheology would not expect to derive their knowledge of each other’s work 
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from this volume. It follows that their accounts of the mathematical 
basis should attempt to be of value to the experimentalist. "Two deductions 
follow, the first being that, wherever possible, the end result of a theoretical 
calculation should be in a form amenable to experimental verification or 
use, with any limitations imposed by the theory clearly indicated. ‘The 
second is that the presentation should be one which does not introduce, for 
reasons of brevity or mathematical elegance, notations or steps which are 
needlessly difficult to follow, and which are not of general value. This 
does not mean that the real difficulties of some parts of the mathematical 
theory can be avoided, but only that it should not be written in a form 
suitable only for fellow-mathematicians. 

R. S. Rivlin, in his chapter ‘‘ Large elastic deformations ”’, has given 
an object lesson in how to present the theory satisfactorily. He has used a 
notation free of tensors, equations being written out in full in terms of the 
co-ordinates. Although this renders them bulky and rather forbidding, a 
little patience will enable almost any research worker to follow his argument. 
He solves particular problems and even indicates where the experimental 
checks of the theoretical results may be found. In contrast, the chapter by 
J. Riseman and J. G. Kirkwood, ‘“ The statistical mechanical theory of 
irreversible processes in solutions of macromolecules ”’, would defeat most 
experimentalists, and might need rather careful study even by fellow 
applied mathematicians. It is therefore out of place in the present work, 
however important a contribution it may be to its subject. 

‘The absence of an adequate account of the bases of classical elasticity 
theory and of the hydrodynamic treatment of Newtonian viscous flow is 
presumably deliberate, although it might have been of some help to 
experimentalists if well carried out. 

The remaining chapters, apart from an Introduction by Eirich, are: 
Phenomenological macro-rheology”’ by M. Reiner, “ Finite plastic 
deformation’ by William Prager, “ Stress-strain relations in the plastic 
range of metals—experiments and basic concepts”? by D. C. Drucker, 
‘“ Mechanical properties and imperfections in crystals” by G. J. Dienes, 
‘‘ Dislocations in crystal lattices’ by J. M. Burgers and W. G. Burgers, 
‘Mechanical properties of metals’ by J. Fleeman and G. J. Dienes, 
‘“Some rheological properties under high pressure”? by R. B. Dow, 


oe 


‘Theories of viscosity ’’ by A. Bondi, ‘‘ Dynamics of visoelastic behaviour ”’ 
by ‘Turner Alfrey and E. F. Gurner, “ Viscosity relationships for polymers 
in bulk and concentrated solution ” by ‘I’. G. Fox, Serge Gratch, and S. 
Loshaek, ‘‘ Streaming and stress bi-refringence’”” by A. Peterlin, ‘* Non- 
Newtonian flow of liquids and solids” by J. G. Oldroyd, and “* Acoustics 
and the liquid state’ by R. B. Lindsay. Inevitably these are uneven in 
merit but a generally high standard has been maintained. 

While the volume indicates that the complete work will fall short of 
rendering other texts unnecessary, it will be a valuable stand-by for all 
working on this subject. Unfortunately the cost per volume is so high 
that it will deter many individuals from acquiring the whole set. 


A. G. WARD 
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